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ANALYSIS OF FINITE-SIZE PHASED ARRAYS OF 
CIRCULAR WAVEGUIDE ELEMENTS* 

By M. C. Bailey 
Langley Resejirch Center 

SUMMARY 

A derivation is presented for the calculation of the interelement mutual coupling in 
a finite- size planar array of waveguide-fed apertures covered by a multilayered dielectric 
and/ or plasma. The general mutual admittance expression is evaluated for circular aper- 
tures and the mutual coupling calculations are verified experimentally for two transverse 
electric (TEji) circular waveguide mode excited apertures. A parametric study of higher 
order mode aperture fields indicates that the only significant change in the circular aper- 
ture mutual coupling is due to the transverse magnetic (TMu) mode, which introduces an 
additional phase shift. Qualitative agreement between calculations for a 183-element array 
of circular apertures and an infinite array establishes the validity of the finite- array theo- 
retical model. 


INTRODUCTION 

The wide flexibility available in the design of antenna arrays is very useful in appli- 
cations where factors such as beam shaping, side lobe level control, and rapid beam steering 
are of prime consideration; however, the design is complicated by the effects of mutual 
interaction between the radiating elements. These interactions are principally evident as 
(1) a distortion of the radiation pattern, (2) an element driving impedance which varies as 
the array is phased to point the beam in different directions, and (3) a polarization variation 
with scan angle in an array with elements which can support more than one sense of polari- 
zation. The degree to which the interelement coupling affects the performance of the array 
will depend upon the element type, the polarization and excitation of each element, the geom- 
etry of the array, and the surrounding environment. In order to study the effects of mutual 
interelement coupling in an array, the analysis must include all these factors. 

The work reported here is an analysis of the mutual coupling in a planar array of 
circular waveguide-fed apertures in an infinite conductor as typically illustrated in figure 1. 


*The information presented herein was offered as a thesis entitled "Near Field 
Coupling Between Elements of a Finite Planar Array of Circular Apertures" in partial 
fulfillment of the requirements for the degree of Doctor of Philosophy in Electrical 
Engineering, Virginia Poljdechnic Institute and State University, Blacksburg, Virginia, 
December 1972. 




The analysis is expected to yield good results for planar arrays on finite-size ground 
planes which are electrically large; however, for small ground planes or for array ele- 
ments near the edge of a finite ground plane, the electromagnetic scatter from the ground 
plane edges may be significant in some cases. 

The problem is first formulated for arbitrary waveguide apertures radiating into 
a multilayered region and then specialized to circular apertures excited in either TE or 
TM modes. The effects of mutual coupling are determined by first computing the self and 
mutual admittances among all the elements of the array to form a complex admittance 
matrix, which is then operated on to determine the complex scattering matrix for the 
array. The scattering matrix gives the relationship between the amplitudes and phases of 
the waveguide modal fields which are incident on and reflected from the apertures. This 
relationship then allows one to determine the reflection coefficient and coupling coefficients 
of all the elements of the array for any phasing or amplitude taper. 


A 

A or A(xi,yi,Zi) 

A(k^, ky, Zj) 
A'(k^, ky, Zj) 

A(a, is ) 


A(a,/3, 0) 


SYMBOLS 

magnetic vector potential 

functional form for z^ component of A 

bidimensional Fourier transform of A(x^ y., z^ 

derivative of A(kj^, k^, z^ with respect to z^ 

undetermined quantity used in equation (69) 

bidimensional Fourier transform of A(Xj^, y., z^ in cylindrical 
coordinates evaluated at z^ = 0 
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[a] 

B(a, /3) 


[b] 

C(a,yS), Cj(a,/S), C2(a,/3) 


^TM 

D(a, /3) 
d 

^2> • • •’ ^p-1’ *^pA 
dp+1, • • %'-!> %'j 

E 


parameter defined by equation (79) 

parameter defined by equation (80) 

dummy parameter used in equations (132) and (133) 

radius of ith circular aperture 

radius of jth circular aperture 

complex amplitude of pth waveguide mode incident on ith 
aperture 

complex amplitude of qth waveguide mode incident on jth 
aperture 

complex column matrix whose elements consist of all ap 

undetermined quantity used in equation (69) 

complex amplitude of pth waveguide mode reflected from ith 
aperture 

complex column matrix whose elements consist of all the bp 

undetermined quantities used in equations (70), (53), and (54), 
respectively 

quantity defined by equation (81) 
quantity defined by equation (82) 
undetermined quantity used in equation (70) 
thickness of one dielectric layer 

distances from aperture plane to outer surfaces of layers 
1, 2, . . ., p-1, p, p+1, N’-l, N’, respectively 

electric field vector 



E(xj, yj, 
E(k^, ky, Zj) 


F 

F or F(Xp y., 
F(k^, ky, Zj) 
F'(k^, ky. z.) 

£(a, / 3 , Zj) 
f'(a,/3, Z.) 
g(a, /3, Z.) 
g'(a, / 3 , Z.) 

H 

H(x., yp z.) 

H(k^, ky, Zj) 


functional form of E 

bidimensional Fourier transform of E(Xj^, y^, z^ 

normalized electric vector mode function for qth TE waveguide 
mode 

normalized electric vector mode function for qth TM wave- 
guide mode 

electric vector potential 

functional form for component of F 

bidimensional Fourier transform of F(x^, y^, z^) 

derivative of F(k^, k^, z^ with respect to z. 

normalized form of F(kj^, k^, z^) defined in equation (56) 

derivative of f(a, /?, z^ with respect to z^ 

normalized form of A(k , k , z.) defined in equation (57) 

X y 1 

derivative of g(a, /3, Zj^) with respect to z^ 
magnetic field vector 
functional form of H 

bidimensional Fourier transform of H(xp y^, z^^) 

normalized magnetic vector mode function for qth TE wave- 
guide mode 

normalized magnetic vector mode function for qth TM wave- 
guide mode 

reaction integral (eq. (21)) 
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Il 3 I 9 J 




J' (z) 
m' ' 


f-r 


k' k' 

Kjj, Ky 


m., m., m!, m! 


intermediate quantities used in derivation of equations (167) 
to (170) 

equivalent current for pth mode in ith aperture 
equivalent current for qth TE waveguide modal fields 
equivalent current for qth TM waveguide modal fields 
complex column vector whose elements consist of all I 

Pi 

Bessel fvinction of the first kind of order m and argument z 
derivative of Jjj^(z) with respect to z 

wave propagation constant in free space, 2v/\ 

Fourier transform variable with respect to 
Fovirier transform variable with respect to 
complex wave propagation constant in z^ direction 
dummy variables for integration 
cutoff wave number for qth TE waveguide mode 
cutoff wave number for qth TM waveguide mode 
number of waveguide modes in each aperture 
total number of modes assumed in the jth aperture 
order of Bessel functions and cyclic variation of fields 
number of apertures in array 
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N' 


R 


Si 


Pi,qj 


[S] 


T T 

^x’ V 


t 

TE 

TM 


Vi^(/3) 


V, 



V" 

q 


number of dielectric layers outside of aperture plane 

center-to- center spacing between two apertures or between 
origins of ith and jth aperture coordinate systems 

area of ith aperture 

complex coupling coefficient between pth mode in ith aperture 
and qth mode in jth aperture 

complex square matrix whose elements consist of all Sp ^ 
beam-pointing directional cosines 
time, sec 

transverse electric 
transverse magnetic 

quantity for simplification of admittance expression (see 
eqs. (167) to (170)) 

quantity for simplification of admittance expression (see 
eq. (167)) 

equivalent voltage of ith aperture field 
equivalent voltage of jth aperture field 
equivalent voltage for pth waveguide mode in ith aperture 
equivalent voltage for qth waveguide mode in jth aperture 
equivalent voltage for qth waveguide mode in kth aperture 
equivalent voltage for qth TE waveguide modal fields 
equivalent voltage for qth TM waveguide modal fields 
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[VI 


w 


X, y, Z 


*i’ ■’!> “i 


yj’ 


X, y, Z 


r •'1’ 1 


yj’ "i 


x|,y! 


^i’^i 




Y 

Pi,qj 

[Y] 

[Yq] 

z’ 


complex column matrix whose elements consist of all V 

Pi 

quantities defined by equations (162) and (163) 
dummy variable 

variables in reference Cartesian coordinate system 
variables in ith aperture Cartesian coordinate system 
variables in jth aperture Cartesian coordinate system 
unit vectors in x, y, and z directions 
unit vectors in Xp y^ and z^ directions 
unit vectors in Xj, y^, and z^ directions 

translation of Xp yp coordinate system in x and y directions 


translation of x^, y^, z^ coordinate system in x and y directions 

characteristic admittance of pth waveguide mode in ith aperture 

element in ith row and jth column of [y] or mutual admittance 
between ith aperture electric field and magnetic field pro- 
duced by jth aperture field for one mode apertures 

mutual admittance between pth waveguide mode electric field 
in ith aperture and magnetic field produced by qth waveguide 
mode in jth aperture 

complex square matrix whose elements consist of all Y 

Pi,qj 

complex diagonal matrix whose nonzero elements consist of 


dummy variable used in definition of delta function (see eq. (28)) 
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<x 


/? 


3(z - Z*) 
e or e (Zj) 

£’ 

^0 

£l(0) 

X. 

^ or yx(Zj) 

^0 

^j(O) 

^P(fi), ^™(5) 


d,y 


P 


angular Fourier transform variable in cylindrical coordinate 
system of ith aperture 

normalized radial Fourier transform variable in cylindrical 
coordinate system of ith aperture 

delta function defined by equation (28) 

permittivity of dielectric region 

permittivity of medium outside of layered region 

permittivity of free space 

permittivity for = 0"^ immediately adjacent to aperture plane 
quantity defined by equation (165) 

quantity defined by equation (165) with i replaced by j 

wavelength in free space 

permeability of dielectric region 

permeability of medium outside of layered region 

permeability of free space 

permeability for z^ = O'*' immediately adjacent to aperture plane 
quantities defined by equations (164) and (166) 

quantities defined by equations (164) and (166) with i replaced 
by j 

dummy variables 

dummy variable for integration in equations (132) and (133) 
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radial variable in ith aperture cylindrical coordinate system 
radial variable in jth aperture cylindrical coordinate system 

sum of terms for all values of q 

sum of terms for all values from j = 1 through j = N 

sum of terms for all values from Qj = 1 through qj = Mj 

magnetic scalar potential 
angle defined by equation (89) 

angular variable in ith aperture cylindrical coordinate system 
angular variable in jth aperture cylindrical coordinate system 

rotation of x^, y^ coordinates with respect to x, y coordinates 
rotation of Xj, yj coordinates with respect to x, y coordinates 
relative polarization angle between i and j aperture fields 
solution to equation (5) 
n!th zero of 
njth zero of 

electrical scalar potential 

temporary variable used in derivation of equations (167) to 
(170) to represent quantity in equation (95) 


solution to equation (6) 



phase shift between array elements for H-plane scan 


phase shift between array elements for E-plane scan 
angular frequency, rad/sec 




3X? 92^ 


Subscripts: 

i ith aperture 

j jth aperture 

ij either the ith row and jth colvunn of matrix or the interaction 

of the jth aperture fields upon the ith aperture fields 

m order of Bessel function (see eq. (116)) 

m^ first subscript of transverse electric waveguide mode in ith 

aperture and order of Bessel function in field equations 

mj first subscript of transverse electric waveguide mode in jth 

aperture and order of Bessel function in field equations 

m! first subscript of transverse magnetic waveguide mode in ith 

aperture and order of Bessel function in field equations 

mj first subscript of transverse magnetic waveguide mode in jth 

aperture and order of Bessel function in field equations 

N' outermost dielectric layer 

N'+l medium outside the layered region 
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*^1 




p 


P+1 

Pi 

q 

t 

X^, yj, Zi 

Superscripts: 

(i) 

0 ) 


P 


second subscript of transverse electric waveguide mode in 
ith apertvu*e 

second subscript of transverse electric waveguide mode in 
jth aperture 

second subscript of transverse magnetic waveguide mode in 
ith aperture 

second subscript of transverse magnetic waveguide mode in 
jth aperture 

pth dielectric layer, except when used in conjunction with e’ 

^ JT 

or to denote pth waveguide mode functions 
adjacent dielectric layer outside pth layer 
pth mode in ith aperture 
qth waveguide mode 
qth mode in jth aperture 
transverse component of field vectors 
components in ith aperture Cartesian coordinate system 
components in jth aperture polar coordinate system 


electric field in ith aperture 

either the electric field in jth aperture or the magnetic field 
produced at ith aperture due to an electric field excited in 
jth aperture 

pth dielectric layer 
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(Pi) 

TE 

TM 

TE, TE 
TM, TM 
TE, TM 

TM, TE 


pth waveguide mode in ith aperture 

transverse electric waveguide mode 

transverse magnetic waveguide mode 

mutual admittance between TE modes in apertures i and j 

mutual admittance between TM modes in apertures i and j 

mutual admittance between TE mode in ith aperture and TM 
mode in jth aperture 

mutual admittance between TM mode in ith aperture and TE 
mode in jth aperture 


REVIEW OF THE LITERATURE 

There are many approaches to the analysis of mutual coupling effects upon the 
performance of phased arrays and each one has its own inherent advantages and disad- 
vantages. It is not the intention of the author to present an exhaustive review of all the 
previous work that has been accomplished in the analysis of phased arraysj however, a 
summary will be given of the more pertinent work of which the author is presently aware, 
and more specifically that which is applicable to planar arrays of apertures. This summary 
is presented in order to acquaint the reader with the scope, depth, and variety of attention 
which phased arrays have received durii^ the past decade. 

The theoretical analyses can generally be divided into two broad categories such as 
infinite arrays and finite arrays. The infinite-array approach is very useful in the analysis 
of the impedance and radiation characteristics of the elements near the center of a very 
large array, but breaks down when applied to the elements near the edge. The finite-array 
approach yields good results for all the elements of the array, but the analysis is more 
complicated, requires more computer time to obtain results, and is generally restricted to 
arrays of no more than about 200 to 300 elements because of the necessity of inverting a 
large matrix or solving a set of simultaneous equations. 

Much effort has been devoted to the analysis of a variety of infinite arrays of 
periodically spaced identical elements (refs. 1 to 57). These have included the more 
common aperture elements such as infinite slots (refs. 23 to 31), rectangular (refs. 32 to 
47), and circular (refs. 48 to 54) as well as the ridged waveguide aperture (ref. 55) and 
multiple frequency interleaved arrays (refs. 56 and 57). Some authors have also considered 
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the effects of dielectric loading such as plugs in the waveguide apertures (refs. 29, 31, 40, 
42, 51, and 54) or dielectric sheets covering the aperture plane (refs. 11, 21, 24, 26, 28 
29, 30, 40, 41, 49, 51, and 54), and the effects of higher order aperture fields (refs. 31, 

38 to 40, 42, 44 to 46, and 50 to 54). 

These analyses have been very useful in the study of certain resonance phenomena 
which have been observed in large phased arrays and array simulators. (See refs. 1, 4, 

11, 38, 39, 43, and 58 to 64.) This resonance is manifested by a null in the array element 
pattern or a large reflection coefficient at specific scan angles closer to broadside than 
the angle at which a grating lobe can occur; thus, the angular scan range of a large phased 
array is limited. This resonance could be considered as the electromagnetic analogy of 
the Woods "anomalies” (ref. 65) for the diffraction of light from optical gratings. This 
resonance in infinite arrays is generally attributed to the excitation of surface waves on 
the periodic structure (refs. 4, 11, 30, 51, 52, 54, 58, 59, and 66) or higher order mode 
aperture fields (refs. 31, 38, 39, 45, 53, 67, and 68). 

Several techniques are available for the elimination of this resonance or of improv- 
ing the wide-angle matching capability of large arrays (ref. 69). These techniques involve 
the use of such things as conducting fences or corrv^ations between the radiating elements 
(refs. 70 to 74), irises in the apertures (refs. 74 to 77), proper design of the dielectric 
loading (refs. 78 to 81), separate matching networks for each element (refs. 82 and 83), 
interconnecting circuits (refs. 84 and 85), selective mode excitation (refs. 86 and 87), or 
possibly disrupting the periodicity of the array (refs. 88 to 92). The wide-angle matching 
is achieved either by a reduction in the interelement mutual coupling or by proper com- 
pensation. In either case, a detailed knowledge of the interelement coupling or terminating 
impedance is required. 

The theoretical analyses for infinite arrays and measurement techniques (refs. 93 
to 97) have proven useful in the study of the radiation and impedance characteristics of the 
"typical" elements of large arrays; however, the "nontypical" elements near the edge or 
the elements of a small array must be analyzed by other means. 

The characteristics of the edge elements in large arrays have been analyzed by 
perturbation (ref. 98) and modifications (refs. 99 and 100) of infinite- array techniques. An 
integral equation method has been used to study the radiation properties of a finite parallel- 
plate waveguide array (refs. 101 to 103). These studies indicate that the impedance and 
radiation properties of the edge elements of an array can be vastly different from those 
near the center. 

Much effort has also been devoted to the determination of the mutual coupling between 
pairs of waveguide apertures. The most comprehensive study of the coupling between 
various antennas was performed by a group at the University of Michigan (ref. 104); however, 
others have also made significant contributions in this area by using a variety of techniques. 
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Graf (ref. 105) investigated the effect of mutual coupling between half-wave slots by using 
the electromagnetic duality of slots and dipoles. Tartakovskiy and Rubinshteyn (ref. 106) 
introduced a numerical method for solving the system of Wiener- Hopf-Fok equations 
which occur in the diffraction at a finite or infinite number of equidistant half-planes and 
applied it to the coupling between two waveguides. Others (refs. 107 to 111) have used 
Keller's geometrical theory of diffraction (ref. 112) to compute the coupling between 
parallel-plate waveguides. Others (refs. 113 to 121 and 123) have used variational technique 
to determine the mutual coupling between rectangular (refs. 114 to 121 and 123), parallel- 
plate (refs. 121 and 122), and annular slot apertures (ref. 113). Some have also considered 
the effects of a dielectric or plasma outside the aperture plane (refs. 118 to 121). Fante 
(ref. 121) used the concept of an impedance sheet to represent the plasma layer under 
certain restrictions. Galejs (ref. 118) approximated the external plasma layers by a large 
dielectric-filled waveguide. Golden and Stewart (refs. 119 and 120) analyzed the coupling 
between rectangular slots under an inhomogeneous plasma by using an integrated electron 
density and a stepped approximation for the plasma profile. Previous work (ref. 122) has 
indicated that stepped plasma profiles can sometimes yield erroneous resonance effects 
which are not present in a practical plasma. Sugio and Makimoto (ref. 123) formulated a 
variational expression for the scattering coefficients of a finite array of rectangular wave- 
guides with dielectric pli^s; however, no results were given. 

The work to be presented in this paper is a variational formulation for the mutual 
admittance of two waveguide apertures which need not be identical in shape nor excitation. 
The formulation is general enough to include the effect of an arbitrary number of dielectric 
and/or plasma layers, each of which may be inhomogeneous; however, no stepped approxi- 
mation to the plasma profile is made nor is an integrated electron density approximation 
used. 

Since no results have been published for finite arrays of circular waveguides, the 
general formulation for mutual admittance is evaluated for circular apertures excited in 
either TE or TM circular waveguide modes and numerical as well as experimental data 
are presented for mutual coupling with either free space or a dielectric sheet outside the 
aperture plane. 

The approach used in the general formulation parallels that for the self admittance 
of one aperture (ref. 124). 


THEORY 
General Theory 

It is assumed that each aperture in the array is fed by a uniform waveguide, the 
cross section of which coincides with the aperture. The electromagnetic fields in the 
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apertures will be represented as the sum of the waveguide modal fields; therefore, the 
total transverse fields in each aperture are given by 


E. = V! e’ + V V"e" 

t ^ / ■ q q 


H. = r li’ + I" fi" 
t Z_, q q / . q q 


( 1 ) 


( 2 ) 


where e^, h^ and e^, h|j represent the normalized vector mode functions for the TE and 
TM modes, respectively, defined so that in Cartesian coordinates 


e’ = - f— X + y\ s 

V 3x 3y / 




h' = z X e' 

q q 


e" = z X /J_x + ^y\ ^ 
^ \3x 3y / ^ 




h" = z X e” 

q q 


(3) 


(4) 


where x, y, and z are unit vectors in the x, y, and z directions, and and 0 are scalar 
functions which satisfy the differential equations 





(5) 





( 6 ) 


subject to the appropriate boundary conditions of the waveguide modal fields. 
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The equivalent modal voltages and currents are defined as 



Ex * e* dx dy 

I 4 








Ej-«'^dxdy 

‘q'J 


where the integrals are taken over the cross section of the waveguide. 

Because of orthogonality properties of the vector mode functions 



(P = q) 


(P 7^ q) 
(P = q) 
(p/q) 


( 7 ) 


( 8 ) 


( 9 ) 


( 10 ) 


( 11 ) 


energy propagates along a uniform waveguide in each mode independently; therefore, for 
computational purposes, each modal field in each aperture of the array is assumed to be 
fed by a separate waveguide which can only be excited by that single mode. This assumption 
corresponds to treatii^ an array of N waveguide-fed apertures as an N times M microwave 
equivalent network, where M is the total number of modes needed in each aperture to 
represent the total field distribution adequately. This assumption restricts the analysis 
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to apertures of relatively simple shapes (such as rectangular, circular, elliptical, etc.) 
for which the corresponding waveguide modal fields can be determined. 

The transverse electric and magnetic fields of the p^th mode can be represented 
either as the superposition of an incident (ap.) and reflected (bp.) wave, or as an 
equivalent voltage (Vp^ and current (Ip^. For TE modes 


























( 12 ) 


(13) 


where Yp. is the characteristic admittance of the p^th mode. The corresponding 
expressions for TM modes are obtained by replacing the primes by double primes 
in equations (12) and (13). 

Because of the coupling or the mutual interaction of the external fields, the 
equivalent aperture voltages and currents will not be independent, but will be related by 
a set of simultaneous equations such as 



N 




M. 


Pi,qj 


(14) 


where N is the total number of apertures in the array and Mj is the total number of 
modes in the jth aperture necessary to represent the aperture field adequately. 

The amplitudes of the incident and reflected modal fields are related by a similar 
set of simultaneous equations such as 



N 


E 


M, 



PpQj Qj 


(15) 
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If each aperture requires M modes to represent the aperture fields adequately, 
then N times M equations such as equations (14) and (15) would be needed to describe 
the couplii^ mechanism of the array. In matrix notation these equations are written as 


[I] = 

[Y] 

[V] 

(16) 

[b] = 

[S] 

[a] 

(17) 


By algebraic manipulation of equations (16) and (17), the wave scattering matrix [S] will 
be related to the aperture admittance matrix [Y] as ” 

[S] = [Yq] - YJ] [Yq] + Y]]"^ (18) 


where [Yq] is a diagonal matrix whose elements are the characteristic admittances of 
the waveguide modes, and [ ] indicates matrix inversion. Thus, the number of 
^ertures N and/or the number of modes M per aperture is limited by the ability of 
the available computer to invert an N x M complex square matrix. 

The Coupling problem then reduces to the determination of the elements of the 
aperture admittance matrix which are the mutual admittances between each aperture modal 
field and all others of the array. 

Mutual Admittance Between Apertures 

General .- In order to compute the coupling between apertures, the components of 
the admittance matrix must be determined. As seen from equation (14), the component 
Yp q (where p^ refers to the pth mode in the ith aperture and qj refers to the qth 

mode in the jth aperture) is the mutual admittance between modes p^ and q^ with all 
other modal voltages set equal to zero; that is, 


Pi,q 




(19) 


with all V = 0 except V . 

Ok Oj 

In order to simplify the subscript notation, and since each modal field will be 
treated as a separate aperture, the notation Y. . will be used to represent the (i,j)th element 
of the (N times M) by (N times M) admittance matrix. 
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A stationary expression for mutual impedance for linear antennas or its dual for 
aperture self admittance can be obtained from the electromagnetic reaction of the assumed 
equivalent electric or magnetic currents (ref. 125, sections 7-9 and 8-12). The mutual 
admittance between two apertures also can be determined from a consideration of 

Y. . = ^ ff X • z. dS. (20) 

VjV. JJ ^ ^ 

' ^ Si 

where V. and V. are the normalized modal voltages (see eqs. (7) and (8)), is the 

1 ] 

assumed electric field of the ith aperture, H'-’' is the magnetic field produced in 
aperture i by an assumed electric field E^^^ in aperture j. The integral in equa- 
tion (20) is taken over the area (Sj) of the ith aperture. 

Borgiotti (ref. 115) used equation (20) to show that the mutual admittance of identical 
apertures radiating into free space can be expressed as the Fourier transform of a function 
which is obtained from the plane- wave spectrum of the field radiated by the aperture. He 
also showed that this formalism can be used to determine the "grating lobe series" for the 
driving-point admittance of an element in an infinite periodic array of identical apertures. 

A more general expression is developed here which is applicable to apertures 
which are not identical in shape or excitation. The mutual admittance expression will also 
include the influence of a planar stratified region outside the aperture plane as indicated 
in figure 2. This expression, which is not presently available in the literature, is then 
used to compute the near-field coupling between circular apertures in a finite planar array. 

Since the tangential component of the assumed aperture field E^^^ is zero over 
the remaining surface of the infinite aperture plane (all other aperture voltages are tem- 
porarily set equal to zero, see eq. (19)), the surface integral in equation (20) can be extended 
to infinity 


Si 

where x^ and y^ are the coordinate variables of the ith aperture. Taking Fourier 
transforms so that 


[eW X • Zj dSj = J J [e^^^ X • Zj dXj dy^ 


( 21 ) 


-00 -00 


E^^^ (kjj, ky) = f f (x., y^ e^ ^ ^ e^ dx^ dy. 


( 22 ) 
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e' m' 


I 


Figure 2.- Cross section of N waveguides radiating into 
N' dielectric layers. 


SG) (k^, kp = [ r (x., y.) 


H 


jk^i jk'yj 


e " * e ^ ^ dx. dy. 

1 •'i 


-00 ‘'-00 


and inversely, 


( 277 )^ lo, 


s<‘) 


H*'* (xp yi> = J_ r r g‘i’ (k;;, ky 

( 2 v )^ *^-00 *^-00 


(J).W .M/^^"ie'^^idk-dk; 


(23) 


(24) 


( 25 ) 


and substituting equations (24) and (25) into equation (21) gives 
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00 00 00 00 fr /.\ / \ 1 

— I I Jf I I f OSj, ky)xH^j) VJ 

2 77) -00 •'-00 \ y~CO •'-00 *^-00 •'-00 


iVl i’iVi /‘Vi 

o •> Q e •' 


dk^ dky dk^ dkj^ > dx. dy . 


Interchanging the order of integration yields 


i? £ £ I 


(k^, ky) X j I 50) (k^, kp 

rr\ 


1 f“ f” "j^+^pyi I 1-1 

1 J.V >' “iK^ ict 


and using the definitioii of the delta function s(z - z'), (eq. (C-19), ref. 125) 


8(z - z’) = — f dw 

2rr 


equation (27) can be written as 


00 a 

— f f 

(2t7)^ ‘'-CO ■’-a: 


(k , k ) 

V x> y ' 


00 a 

/ / 

•' 00 *'_rr 


(k;^, kp s (k^ + k;^) .s':(ky + kp| dk^ dk^ 


which yields 


•Zi^dk^dk 


(kx> ky) 


X (- k^, - ky 


)] • Zi dk^ 


Equation (30) is recognized as a form of Parseval's theorem (ref. 125, eq. (C-15)). If k 
and k , are the wave propagation numbers in the x. and y. directions, then one could 

y 1 1 

visualize H'-’'' (- k , - k ) as the bidimensional Foimier transform of a wave whose direction 
X y 
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of propagation in the x^, plane is reversed. The problem now reduces to the determi- 
nation of the tangential component of (- 1^, - k^) at the ith aperture because of an 
assumed electric field in the jth aperture. 

The electric and magnetic fields external to the aperture plane can be uniquely 
determined from a set of vector potentials 


A = A(x., y., z.) z. 
F = F(x., y., z.) z. 
as follows (see appendix A) 


r. / \ 1 ^2 /A\ 1 3F 

(Xp y=, zj = 5 _ 

i ]coe Bx.Bz^ \iJ.J € dy^ 


^ ^ 1 ^ /A\ ^ 1 BF 

E„ (x., y., zj = + 

^i ]oje By^Bz. \m/ e Bx. 


1 

3 

"1 

3 

/A\ 

j" 


e 

^i 

\M/ 


jwA 


TT / ^ 1 3 /F\ ^1 3A 

(x., y., zJ = + 

^i ^ ^ ^ jcu/Lt 3x.3z. Ve/ M Byj 


H (xp yp Zj) 


1 3 /F\ 1 3A 


jco/x 3y. 3Zj \e / n 3x- 


Hz (Xj, y^, z^ = J: -i- 
^i ^ ^ ^ i^.3z. 


'1 3 /F 
/i 3Z. \€ /J 


- jci)F 


(31) 

(32) 

(33) 

(34) 

(35) 

(36) 

(37) 


where £ and /x are the permittivity and permeability of the external medium and oj 
is the angular frequency of the sign^ A time harmonic variation of the form e^"* has 
been suppressed. 
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Substituting the inverse Fourier transforms (eqs. (A27), (A28), (A29), and (A30)) 
into equations (32), (33), (35), and (36) gives, after interchanging orders of integration and 
differentiation. 




*V- A'(k^, ky. Zj) - Jk^F(k^, ky, z,) 


«X *X> ‘'y = - iky AOSc, k Zj) F-(k^, k Z.) 




\K’ ky. = l>Sc AOS;, ky, zp - ^ F'(k^, ky, Zj) 


where the primes denote differentiation with respect to z^. 
Then equation (30) becomes 


(38) 


(39) 


(40) 


(41) 


I =. 


j Jk? +k^ 


(27ry 


I I W 

-00 - 00 V ' ' 


where 


+k^ 
_x y 

CO/J.(0) 


F(k^, k, . 0) F’(- k^, - k . 0)^ dk^ dk. 


X’ y 


X’ y’ 


X y 


A'(kj^, ky, 0) = 


_ A(k^, ky, zp 


z.=0 


(42) 


(43) 


F'(- k^, - k 0) = 


— F(- k^, - ky, zp 


z.=0 


(44) 


If all apertures except the jth are short circuited, then continuity of tangential elec- 
tric fields over the aperture plane gives, from equations (38) and (39), 
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A'fk k 0) = - 
A (k^, ky, u; 


F(kj^, k , 0) = 


[■ScE^K- V “> * Vyi V>> 


a) 


j(k^k2) 


^y’ 


(j) 


(45) 


(46) 


where (k , k , 0) and (k , k , 0) are the bidimensional Fourier transforms 


X. ' X’ y 


y. ' X’ y' 

of the assumed modal electric field in the jth aperture. 

Likewise, if all apertures except the ith are short circuited, 


A’(- - k 0) = 


we 


(0) 


2 . ,,2 
y 


k^.k: 


-1 


r (- k ,, - ky , 0 ) = 




Vx^- >Sc> - V - V® <- *Sc’ - V ®>] 


(47) 


(48) 


3(k^-9 

where E^ (- k^, - k„, 0) and E^j^^ (- k^, - ky, 0) are the bidimensional Fourier trans- 


"X. ' 'X' y 


forms of the assumed modal electric fields in the ith aperture with the direction of propa- 
gation in the x. and y^ directions being reversed. 

Note that the transformed wave equations (eqs. (A3 5) and (A3 6)) are even functions 

of kjj and ky; therefore, by using equations (45), (46), (47), and (48) in equation (42), 


the mutual admittance becomes 


Y - 1 r r f 

-jo;e(0) 

"A(kj^, ky, 0)”j f 

VjVj(27r)^ -^-coi 




k (- k , - k , 0) 

y y _ \ X’ y’ 


MO) 04 ■+ ky)j 


F’(k^,ky, 0)^ 




Vxi V Vyj^c'V 


• s^% ‘Sc’ - N’ °) - *Sc4i ‘Sc’ - ‘V’ 


dlScdky 


(49) 
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If a change of variables is made in the transform domain to cylindrical coordinates 
so that = Kq/Scos a and = kQ/3 sin a , 'then 




Yy = — : — ^ I I 

VjVj(277)^ /3=0 0=0 


jA'(a,/3,0) M *i 


; (a, P, 0) COS d + EJ^/ (a, P, 0) Sin 


X Te^^) (a,-P, 0) COS a + E^‘^ (a, - P , 0) Sin a + 

L*i -I ko^F(a,/3,0) 

_ ^0 


E^^ (a, p, 

_ i 


0) sin a - (a, /3, 0) cos a 


X (a, -/ 3 , 0 ) Sina . (a, ♦ / 3 , 0 ) COS a > /3 d /3 da 


where A{a^/3 ^ 0) and F(a,y5 , 0) now satisfy the differential equations 


A/ « i 1 d 

A(a, /3, z.) - _ 


e(z.) dZj dz. 


A O MzJ « (zj 9 

J- A(a, yS, z.) + — i L-/3^ A(a,yS,Zj)-0 (51) 


1 d^(z.) j 9 

^ F(a,^, Z.)-— L _1 F(x,/3, Z ) +ky 

, 2 ^ /x(z.) dz, dz. 1 U 

dz^ ^^1 1 1 


9 ^ ( 2 .:) 9 

0 ^ F(a,^,z^=0 (52) 

^0 ^0 


subject to the boundary conditions (eqs. (A37), (A38), (A39), and (A40)) at each boundary 
(zi = dp) in figure 2. 

Assume that the region outside the aperture plane {z^ > 0) consists of N' layers 
whose total thickness is dj^,. Also assume that the remaining space outside the layered 
region (z. > dj^,) is filled with a homogeneous material whose permittivity and perme- 
ability are €* and m', respectively. The solutions to equations (51) and (52) outside the 
layered region (z. > dj^,) will then be of the form 

“ jk z • 

(a,/3,Zj)=Ci(a,^) e (53) 


(tt > /5 > Zj) — C^{o.y p) e 




where k is defined to satisfy the radiation condition at infinity, that is, 
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kz = - i^o 






e'jjLt 





( 55 ) 


For convenience, the solutions to equations (51) and (52) for each layer p will be 
normalized to the solutions in the outer region evaluated at the outer surface of the layered 
region (z^ = d^^,) according to 


fp(a, - 

gp(a,^, z^ = 


Fp(a,/3, z^ 
^N'+l 

Ap(a, p, 

•^N' +1 (^ » ^ ’ 


(56) 


(57) 


which are solutions of 


d^ . , 1 de(zj dg (a,^,z.) 2 

— go(<^, /3, z.) - L . P 

dz2 £(zp dz. 


dz, 




m(zJ £ (Z.) 2 

1 L - /3 


^"0 


gp(^,/5,zp=0 (58) 


4 

dzf ^ dzj ° 


M(Zj) g (zp ^ ^2 


^o"o 


fp(“,^,Zi) = 0 


(59) 


Then by using the boundary conditions at each interface (z. = d ) 

' 1 p' 


fp(a> fi, dp) - fp^i(a, P, dp) 

(60) 

gp(a, / 3 , dp) = gp^j(a, P, dp) 

(61) 

M„(d ) 
^P + 1 ' p' 

(62) 
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( 63 ) 


e (d ) 

«p'“’ V ' 7 ^) V 

P+1 ^ p' 


starting with the initial conditions 






^N’ 


‘N’ 




/^j ~ ” ^^2 




(64) 

(65) 


( 66 ) 


(67) 


and solving the differential equations (58) and (59) for each layer in turn beginning with 
the outermost layer and working back toward the aperture pl^e (z. = 0), the mutual 
admittance for two assumed aperture field distributions (E^^^ and radiating into 

a plane multilayered region can be determined by performing the following integrations: 

ei(0) ~ 

ko gj(a, /S, 0) 

jgl(“,/3, 0) 



X 


r(i)( 


E^'(a, yS, 0) COS a+ E^^^a, j3, 0) sin 
^i ^i 


E[^^(a, ~/3, 0) cos a + E)^^^(a, -j3, 0) sin 


(i)/ 


^i 


a 


jq(a, yS, 0) 


L '' ^0 


k„ fi (a, y6, 0) 

U A 


E^^^a, yS, 0) sin a - E^J^(a, y3, 0) COS 
_ X. y. 


^ii), 


“] 


X 1 E^%, - yS , 0) sin a - E[^\a, - yS , 0) cos 


L i 


yS dyS da 


( 68 ) 
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Only a limited number of dielectric profiles have so far been investigated whereby 
the solutions to equations (58) and (59) can be expressed in terms of well-known fvinctions 
A few of these solutions are found in reference 126, No attempt is made here to cover 
this class of problems. It suffices to point out that once the available solutions are 
evaluated in the aperture plane (z. = 0), the mutual admittance can then be determined. 

For the most general case, the differential equations (58) and (59) must be solved 
numerically, but for the special case of a homogeneous dielectric layer (s (z.) = e 

Mp(Zi)=AiQ)> the solutions to equations (58) and (59) take the form ^ ^ ^ 


-jkPz. 

3 ^' + B(a,/3) 


.jkPz. 

gp(a, 13, Z^ = C(a, (3) e ^ ^ 




where the unknown coefficients are determined from the boundary conditions (eqs. (60) to 
(67)). If the external region consists of only one homogeneous layer of thickness d, the 
ratios of the functions in equation (68) become 




where e j and e' are the permittivities of the dielectric layer and the medium outside 
the layer, respectively. 
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If the thickness of the dielectric layer is allowed to go to zero, equations (71) and 
(72) reduce to 


^l(O) 

^0“^ Si (a, /S, 0) 

igi(a,/3, 0) 

jf^(a, /3, 0) 

kn “I— 0) 

U Mq i 




(73) 


(74) 


for a homogeneous half space over the apertures. 

If the permittivity e' is real, the radical ~ represents a branch point at 

y "0 

2 e * 

/3 = — ; therefore, to account for this properly in the integration, the radical must be 
^0 

/ 2 €' 2 e' 

/3 - — for > — , which corresponds to the radiation condition 
^0 0 

(eq. (55)). 

If both ej and e' are real (lossless dielectric layer), it can be seen from equa- 
tions (71) and (72) that the integrand in equation (68) will be infinite for discrete values of 
These poles on the real axis of a complex /S -plane correspond to the excitation of 
surface wave modes and must be properly accounted for by residues for the integration on 
/3 in the vicinity of these poles; however, this problem can be circumvented by assuming 
the dielectric to be slightly lossy, and this assumption causes the poles to move off the 
real /S-axis. In most cases, a dielectric loss tangent of 0.001 is sufficient to eliminate 
the numerical integration difficulty near these poles while maintaining a 3 or 4 significant 
figure accuracy when compared with calculations for a lossless dielectric. 

Circular apertures .- If the apertures are round holes, the fields in the apertures 
can be described by the set of circular waveguide modes whose transverse electric fields 
(normalized according to eqs. (7) to (11)) are given (ref. 127) 
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For TE: 


For TM: 


where 


m.J (A.p.) 

Sin(m-) 

J ] ] A./). J ^ 


(75) 




(76) 


g(j)TM ^ _ ^TM 


pj ' ]' y J m: ' yy 


(77) 


*i> ' '=i 


m! J , (Alp.) 
TM J “j J ^ 


sin (m! d) 
Aj.j J ) 


(78) 


A. = 




(79) 


^!n! 
A* - J J 
] — 


(80) 


cp 


''i 


TE / 






mjnj) 


v2 2 

) -m. 


(81) 


^TM yj^ 

^TM _ j V TT 


^ ^j'^mj+l (x^,„.) 


(82) 
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with 


e =1 

m. 


(m. = 0) 


e = 2 


(m. 0) 


e , = 1 


(m! = 0) 


e 1=2 

“i 


(m! ^ 0) 


and where are the zeros of the Bessel function of the first kind, that is, 


•^m! ^^mjnj) ~ ° 


and y* „ are the zeros of the derivative of J (y) 
''m-n. 

J J J 


J' (x) . =0 

i ^ m.n. 

•’ 1 1 


where the prime on Jjjj (x) indicates differentiation with respect to the argument. 
From figure 3, a transformation of variables is made so that 


X. = R cos (p + cos (</>. + (^ ) 

^ J J 

y. = R sin 4> + p . sin ((^. + 4>^ 


where 


R = y (yj - yp^ + (xj - x!)^ 
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Figure 5*- Coordinate geometry for the ith and jth elements 
of a planar array of circular waveguide-fed apertures. 


is the center-to-center spacing between the apertures, 


^p=^j-^i (88) 

is the polarization angle of the fields in the jth aperture relative to those in the ith aperture, 
the angle 4 > is defined as 


4 > = arc tan 


and the aisles 4 >'^ and </>! are the polarization aisles of the ith and jth aperture fields 
with respect to a fixed x, y coordinate system. 

i 

j 

I 
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are 


From figure 3, the and components of the aperture fields in the jth aperture 


^ V “ ^ j’ ■*■ V 

(Pj, <^.) sin i 4 >^ + 0 p) + E§^ Co ., 0 .) cos (<^ + <^.p) 

Then with the change of variables in equations (85) and (86), and the definitions 
kx = kQ/fl cos a and ky = kp /3 sin a , 


fSc^i ^y^i] " \ ^0^^ cos (a - 0) + kQ/3p. cos 0. - (a - 


and the transforms of the aperture fields (eq. (22)) become 


i/3, a) = 

^i 


a. 2 tt 

• 3 r 


= e^^ r r 

_n 


{.j=o -*5=0 




1 jkn/3p. COS r0.-(a-</> )]■] 

E,^^^ Coj,0j)sin(0j+0p) e^ ^ j 


a. 2 tt 

r f 




(A) 1 jko/?P; COS nfc.-(a-0 ; 

+ Co., <^.) COS (^>j +0p) e ^ ^ L 1 p. 




where 


i//= kg/SR cos {a -4>) 
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Then for TE 


modes, 


( a. rm.J^ (A.p.) 2tt 

Jf’ '“i cos (. 
‘ Vo L -^0 


jkQ/3p. COs|V>.-(a-(^p)J 


pjdp. 


f 


sin (<^. +<^p) cos (m.<?!>j) 


JkQ/3p. COs[0.-(a-<^p)] 


d0. p.dpj 


(^, a) = cp e^^Jj 


a. [m. J„ (A. pj 2 t7 
r 3 3 3 3 r 

J sin (<^. + A) sin (m. 0.) 

'o L A,p, J„ J P > > 


jkQ^pXOS - ( a- 0p)J 


d0. p.dp. 


r^j r ^ 

J '^m. j V 


JkQ^Sp. COS|^(^.-(a-<^p)] 


d^j Pjdpj 
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and for 


™m!n! 
3 3 


^(j)TM 


(^,a) = 


-C™ eW 




“Fmj (A]'’)) -2. 




I Sin (</>. + <^p) sin (mjt^.) 


*^0 


X e 


jkg/Spj cos [<^.-(oL-<f.p)] 


d^. 


^3 


f 

•Vi 


2tt 


3 *Y) 




<^4>i 


Pi 

] ] 


(98) 


E«>™ (A«) = C™ eW 


-“i 




■mjJm, (AjPj) > 


AjPj 


I COS (< 5 !.. + <^.p) sin (m! 0 .) 


10 


X e 


JkQ/fipj COs[0j-(a-</)p)J 


d0. 


p. dp. 

3 3 


^3 

- [ fAjfj) J ®‘" (*i * V <“i*i> 

*'n .1 ‘'a 


jkQ/Sp. cos [</>.- (a-c^p)] 

X 0 


d#J 


P] dpj 


(99) 
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Then by using trigonometric identities, the integrals over in equations (96), (97) 
and (99) can be expressed in terms of integrals of the form 


/ 


I, = I sin 4>^ sin mc^. e 


jkQ/3p. cos 


d</>. 


277 

= f Si. 


I2 = I sin 4>. cos m<3f)j e 


jkQ/3p. COs|^.-(a-,^p^ 


d<?b. 




I„ = cos <f). sin m<^. e 
o J J 


jkfl cos 


<J0j 




= I cos cos m0j e 


ikp 0p. cos [♦j-(«-*p)] 


d^j 


where m = m. for TE modes or m = m! for TM modes. Then 
1 1 


EO)TEfe.)=cTEei* 

^i ^ 4. 


mjJ„ (AjPj) 

J . (I3 COS (^p - sin <^p) 




®2 ?°= *p+l4®‘" V 




.0)TE,^_.),cTEei* f’ 




A,Pj 


(Ij COS <jbp + 13 sin <^p) 






(98), 


( 100 ) 


( 101 ) 


(102) 


(103) 


(104) 


(105) 
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e())TM , 


-C™ ei* p 


L-^ 


(Ij COS <i>p + 13 sin (^p) 






(106) 


E<i>TM = C™ e>* h 


m;Jm.(Aj^j) 

_-^p— 


(I3 COS <^p - sin 0p) 




p. dp . 

] J 


(107) 


By writing the trigonometric functions in equations (100) to (103) as exponentials, 
the integrals on 4>^ become 


ii = -ia5-%-'7*W 


(108) 


^2 ■ ^ ^5 ^6 ■ ^7 ■ ^ 8 ^ 


(109) 




(110) 


^4 “ ^ ^^5 + ^6 ^ ^7 ^ ^ 8 ^ 


( 111 ) 


where 


j(m+l)(a-,^p) 277-(a-^p) jkQ/3p. COS0 

I5 - e I e e 

-(a-(^p) 


d0 


( 112 ) 




de 


(a-'^>p) 


(113) 
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J / . ^ \ 


d& 


(114) 


•^-(a-V 




d^ 


(115) 


(-^p) 


where a change of variables has been made so that & = - (a - <Pp). The integrals in 

equations (112) to (115) are recognized as a form of the Bessel function of the first kind, 
that is (see p. 367, ref. 128) 


2n- y 

.-m n 

J 

-y 


Jm5 gjz cos B 


(116) 


where y is any arbitrary angle. 

Then with the relationship (p. 128 of ref. 129) 


J.m(-) = (-1)“ Jm(^) 


(117) 


equations (112) to (115) become 


m +1 j(m+l)(a- 0 j 


(118) 


1m ml - j(m-l)(a-0 ) 

Ig = 2^(j)^-“ (-l)""'^ ® 


(119) 


m-1 , . j(m-l)(a-<^J 

17 = 2770)”" ® ^ 


( 120 ) 


m 1 m4-1 -j(m+l)(a-<^ 1 

% = 2.0)-”-* (-1)”** J„^l (k„/3,>j) e p' 


(121) 


Substituting equations (118) to (121) into equations (108) to (111) and combining terms 
yields 
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Ij = l^m+1 ^^0 + !)(“- ■5^p)] 

+ '^m-1 j) - 1) (<1 - ^p)]^ 

I 2 = 77(j)“'^^ l^m +1 |(m + 1 ) (a - <^p)] 

'^m -1 P”^ - 1 ) (^ - <?^p)]J 


I3 - I Jj^+1 (kp ^Pj) sin pm + 1 ) (a - <^p)J 


■*■ '^m-1 pm - 1) (a - 0p)jj 

I 4 = 77 ( 3 )“+^ l^m+1 [(m + 1 ) (a - 0 p)] 

“ '^m -1 (^0 pm - 1 ) (a - c^ppj 


By substitutii^ equations (122) to (125) into equations (104) to (107) and using the recurrence 
equations for Bessel functions 


J.v,_i(z) = 


mJ,v.(z) 


+ J’ (z) 


(z) 


mJ^(z) 


- Jm(^) 


the transforms of the aperture electric fields become 
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E«)TE (ft ,) = cp el* 


X < sin 


in T(m. + 1 ) a - m. 0 p 1 J J^.+i(AjPj) j Ocq/^p.) p- dp. 
^ -* •'n 3 J 


•'0 

a. 


- Sin [(mj - 1) a - mj^p] J ' J . (AjPj) , (k^/iPj) Pj dpj 

J J 


E®'^^ 03 , C^ei^ 


+ cos 


r T 

[(m. + 1) a - m. 0 p] J J 1 (A. p.) J 1 (k(j ^p.) p. dp. 

0 J J 

a. ^ 

[(mj-l). -m.*j] J l(AjPj) J i(ko/*V^l 

0 J J 


E«>™Ce,.)=P(j)“i*'c™''" 

Xi J 


X <cos 


[(m! + 1 ) a - mj^p] J (Ajp.) J^. ^ (k^ /3p.) p. dp^ 

n J J 


‘'O 

a. 


+ cos 


[(mj - 1) a - mj<Pp] I ’ 1 (Ajpj) (ypj) Pj dpj 

n J J 


m!+l 


E®™(/ 3 ,p)=paP C™ei* 

J 

[(mj + 1) - mj0 ] J^,^l(A!p.) (k^/Sp ) p. dp. 

*“ J J 


x<sin 


‘'0 

a. 


- sin [(m! - 1) a - m!0p] f ^ J^,_l(A!p.) J^,_i(k(j^Pj) P^ dp. 

0 ^ ] 


(128) 


(129) 


(130) 


(131) 
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The integrals over p. in equations (128) to (131) can now be evaluated in closed 
form (see p. 146 of ref. 129) 
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and for the TM 


modes 






(138) 
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or 


The transform fields in equations (138) to (141) are even functions of /? when mj 
mj is odd and pdd functions of /3 when or mj is even; therefore, 


Xi X. 

(142) 


(143) 

Xi X. 

(144) 

yi 

(145) 


where (a,/3), (a,;fi), (a,/3), and (a, /3) are obtained from 

Xf yj Xj Yj 

equations (138) to (141) by setting \p =<^p = 0 and replacing the j subscripts by i. This 
will complete the evaluation of the Fourier transforms of the aperture electric fields of 
an open-end circular waveguide. 

Since the ratios of the solutions to the wave equation and their derivatives are 
independent of a (see eqs. (71) to (74)), the mutual admittance between a TEj^^ mode 
in the ith aperture and a TEm^^nij mode in the jth aperture can be expressed as 
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where 
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and by replacing m. by m! and m. by m! in equations (147) to (150), the mutual 
admittance between a mode in the ith aperture and a TM , , mode in the 

jth aperture becomes ^ ^ I j j 
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By replacing by mj in equations (147) to (150), the mutual admittance between a 

TE mode in the ith aperture and a TM . , mode inthejth aperture can be written as 
mm. mm. ^ 
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ij 
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and the mutual admittance between a mode in the ith aperture and a TE^ 

mode in the jth aperture becomes (with m^ in equations (147) to (150) replaced by m!) 
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In the evaluation of the integrals on a , equations (147) to (150) become 
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and for 0, a change of variables is made such that d = a - <p, and equations (147) 
to (150) can be expressed in the form of a Bessel function (see eq. (116)) 
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By substituting equations (154) to (161) into equations (146), (151), (152), and (153), 
and using the following definitions: 
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yTE,TE(^) ^ (_ 1 ) j (kQ^R) cos [(m. +mp - m.<^pj 
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where 
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The integration on /3 in equations (167) to (170) must be numerically evaluated. A 
computer program has been written for the evaluation of these equations for the mutual 
admittance of two circular apertures radiating into a multilayered region of up to four 
layers, two of which may be inhomogeneous normal to the aperture plane. A listing of the 
computer program is included as appendix B. 
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DESCRIPTION OF EXPERIMENT 


Hardware was constructed and an experiment was performed for the verification of 
the theoretical analysis. The verification was accomplished by comparing the measured 
and calculated TEjj mode mutual coupling between two circular waveguide-fed apertvu’es 
for various combinations of frequency, spacing, and polarization. The hardware which 
was constructed and assembled for this purpose is shown in figure 4. 



L-72-7617 

Figwe Jj-.- Experimental model for mutual coupling measurements. 

The hardware in figure 4 consists of a 30.48 cm by 60.96 cm (12 in. by 24 in.) flat 
aluminum plate with two 3.81-cm-diameter (1.5-in.) circular waveguide-fed holes which are 
equally distant from the center of the rectangular plate. The circular waveguide sections 
are connected to standard coaxial adapters (RG 50/U rectangular to type N coaxial) by 
25.4 cm (10 in.) circular to rectar^ular linearly tapered transitions. One of these transi- 
tions had been used in a previous experiment (ref. 130) and performed satisfactorily over 
the frequency range of interest. Since the other transition is dimensionally identical, it can 
be expected to give a similar performance. 
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Swivel flanges are used to connect the circular waveguide sections to the tapered 
transitions. This connection allows the polarization of each apertxu'e to be changed by 
rotating the adapter and transition through the desired angle. The electric-field polari- 
zation of both apertures in figure 4 is vertical, as indicated by the position of the excita- 
tion probes of the coaxial to waveguide adapters. 

The circular waveguide sections are flange mounted to the aluminxim plate as 
illustrated in figure 5 by the unassembled cross-sectional view of the aluminum plate and 
one waveguide. The back side of the aluminum plate is recessed and the waveguide end 
extends out past the flange an equal amount in order to maintain accurate alinement of 
the waveguide with the circular hole. By mounting each waveguide to the aluminum plate 
in this manner, the same waveguide assembly can be used with a variety of flat plates with 
different hole spacings. The one shown in figure 4 is for a center-to-center spacing of 
6.35 cm (2.5 in.). Other plates were constructed with center-to-center hole spacings of 
8.89, 12.70, and 17.78 cm (3.5, 5.0, and 7.0 in.) but are not shown since they are identical 
to the one in figure 4 except for the different hole separations. 



Figure 5>- Cross section of experimental model 
illustrating method of mounting circular 
waveguide to aluminum plate. 


All parts for the experiment were purchased as commercial stock items, with the 
exception of the rectangular plates which were machined from stock aluminum. 

The mutual coupling was measured by exciting one waveguide at the coaxial adapter 
and comparii^ the received signal level at the other coaxial adapter to a known reference. 
This measurement was accomplished by connecting coaxial cables to a signal generator and 
a receiver. The other ends of the cables (shown in fig. 4) were then connected together to 
obtain a reference signal level at the receiver. The mutual coupling was then measured 
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by connecting the cable ends to the coaxial adapters and adjusting in-line calibrated attenu- 
ators until the received signal level was the same as the reference level. The net change 
in the calibrated attenuators is observed as the relative power coupled from the input 
terminal at one coaxial adapter to the output terminal at the other coaxial adapter through 
the open ends of the circular waveguides mounted to the flat aluminum plate. 

The insertion loss due to the waveguide assembly was measured by the same method 
with the flat plate removed and the ends of the circular waveguides connected together as 
shown in figure 6. The insertion loss was measured at each frequency of interest and all 
the measured data presented in the section "Results and Discussion" have been corrected 
for the waveguide assembly insertion loss at each measurement frequency. The insertion 
loss correction was between 0.1 dB and 0.5 dB over the frequency range. 






L-72-7616 

Figure 6.- Experimental model for determination of waveguide assembly 

insertion loss. 

RESULTS AND DISCUSSION 
Comparison Between Measurements and Calculations 

The data in this section are presented primarily for verification of the theoretical 
analysis. Data for various combinations of aperture spacings and polarizations are pre- 
sented as a function of frequency in order to test all facets of the analysis. 

The calculations in this and the following sections were obtained from the computer 
program listed in appendix B. The calculated mutual coupling values were taken from the 
appropriate off-diagonal terms of the scattering matrix. 
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Figures 7 and 8 show a comparison between calculations and measurements for the 
four different aperture separations. There is excellent agreement between the measured 
and calculated values in figure 7 for coupling in the E-plane; however, larger variations 
in the measured data as a fvmction of frequency were observed for the H-plane coupling 
in figure 8. This oscillatory variation with frequency is typical of impedance measure- 
ments for antennas with a truncated ground plane. A slight variation can also be observed 
in figure 7; however, it is more pronounced in figure 8 because of a combination of two 
things. First, when the polarization is such that coupling occurs in the H-plane (fig. 8), 
the E-plane dimension of the ground plane is only 30.48 cm (12 in.) as compared with 
60.96 cm (24 in.) for the E-plane coupling of figure 7. Past experience has shown that the 
dimension of the ground plane in the direction of the aperture electric-field vector has the 
most predominant diffraction effect. Second, the mutual coupling in the H-plane is naturally 
lower than in the E-plane and therefore the measurements become more sensitive to 
scattering from objects such as the edges of the ground plane. This later reasoning is 
verified by the data in figure 8 for the 6.35- and 8.89-cm (2.5- and 3.5-in.) spacings for 
which the mutual coupling is stronger and the scatter in the measured data is less. Since 
the 6.35-cm (2.5-in.) spacing yielded good agreement for both polarizations, the remaining 
data in this section are restricted to this spacing. 



Figure 7*“ mode mutual coupling between two circular 

waveguides radiating into free space (E -plane coupling) . 
Dimensions are in centimeters (inches). 
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Figure 8.- TEj_2 mode mutual coupling between two circiilar 
waveguides radiating into free space (H-plane coupling). 

Dimensions are in centimeters (inches). 

Figures 9 and 10 demonstrate the validity of the theoretical analysis for an arbitrary 
polarization of one aperture field with respect to the other. There is some scatter in the 
measured data due to ground plane edge effects; but, in general, the agreement is very good. 

One thing to be observed by the data in figures 9 and 10 is that both the measured and 
calculated results indicate a trend toward complete isolation (Mutual coupling = -oo dB) 
for orthogonal polarization = 90°); however, this is not always true, as shown in 
figure 11. Here the principal electric fields are orthogonally polarized; however, both 
the measured and calculated results show an appreciable level of coupling between the 
aperture fields. These results indicate that, for certain geometries, the cross-polarized 
fields may have a significant influence upon the performance of a large array. Earlier 
analyses of infinite arrays (refs. 51 and 53) have shown some appreciable charges in the 
radiation characteristics of a phased array due to cross-polarized fields when the array 
beam is scanned at an angle far off the axis. 
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It should be noted that the excellent agreement between the measured and calculated 
results in figure 11 demonstrates the feasibility of using the present analysis to study 
cross-polarized effects in finite-size circular waveguide phased arrays. Such informa- 
tion is important in the design of arrays of circular polarized elements, in which case 
the axial ratio or polarization may vary drastically as a function of scan (ref. 51). 

The data in figures 12 to 15 are presented to justify the theoretical analysis for 
dielectric-covered circular apertvires. The dielectric constant (2.6) and loss tangent 
(0.006) used in the calculations are those measured by Von Hippel (ref. 131). The dielectric 
sheets used for the measurement were the same size as the ground plane and the thick- 
nesses were such that only one surface wave mode can exist (ref. 132). 



Figure 11.- TE^i mode free space mutual coupling between 
orthogonally polarized aperture fields. Linear dimen- 
sions are in centimeters (inches). 
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Very good agreement between measured and calculated values was obtained for 
coupling in the E-plane (figs. 12 and 14); however, large variations with frequency 
occurred in the measured data for coupling in the H-plane (figs. 13 and 15) because of the 
reflections of the surface wave from the ends of the dielectric sheet. Previous work 
(ref. 133) has shown that for thicknesses such that only one surface wave mode exists, 
the E -plane dimension of the finite dielectric sheet produces the largest perturbation in 
the measured data. This perturbation is also evident in the present data as observed in 
figures 12 and 14 where, as a result of the larger E-plane dimension, the oscillations are 
smaller and more closely spaced. 



Figure l4.- TE 3_2 mode mutual coupling between two dielectric- 
covered circular wavegtiides (dielectric constant, 2.6,- 
loss tangent, 0.006^ E-plane coupling). Linear dimensions 
are in centimeters (inches). 
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Figure I5.- TE3_]_ mode mutual coupling between two dielectric-covered 
circular waveguides (dielectric constant, 2.6; loss tangent, O.OO6; 

H-plane coupling). Linear dimensions are in centimeters (inches). 

Although large variations in the measured data were observed in some instances, the 
primary differences between the measured and calculated results are due to the finiteness 
of the ground plane and the dielectric sheet, which was not considered in the theoretical 
model. However, the overall qualitative comparison between the measured and computed 
results has established the validity of the theoretical analysis. 

Computed Higher Order Mode Effect 

A brief parametric study was performed in order to determine the influence upon the 
mutual coupling due to higher order modes in the apertures. These data are summarized 
in table I. 

The mutual coupling was first computed for two 0.7 5 -wavelength -diameter circular 
waveguides in which the aperture field distributions were assumed to be that of the TEn 
mode. Next the 4 by 4 complex scattering matrix was computed for two waveguides with 
two modes in each (TEu plus one higher order mode). Then the appropriate value of the 
scattering matrix for the coupling between the TEn modes (S 13 in this case, was com- 
pared with that computed previously for only the TEn mode assumption. These differ- 
ences for both amplitude and phase are listed in table I for several of the next higher 
order modes. 
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TABLE L- CHANGE IN FREE-SPACE MUTUAL COUPLING BETWEEN TExi MODES 
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It is obvious that those modes, whose first index numbers are different, do not 
influence the mutual coupling calculations. Also, the only mode which has any noticeable 
effect is the TMjj and then primarily only in the phase. In a phase scanning array, 
this chaise in phase due to the presence of the TMjj mode may be of some significance 
and probably should receive more attention in future work. 

The calculations, obtained by assuming that the aperture distributions contain the 
first four modes (Sjg of the 8 by 8 matrix) whose first indices are the same, are also 
compared with the calculations obtained by assuming only the mode. It should be 

noted that the inclusion of additional modes other than the TEjj and TM^j has a 
negligible influence upon the mutual coupling. 

Phased-Array Calculations 

Data are presented in this section to illustrate the variation of the reflection coef- 
ficient as a function of scan for the elements of a finite planar array of circular aper- 
tures excited in the TEii mode. The elements of the array will be in an equilateral 
triangular grid arrangement as indicated in figure 16. The dimensions were chosen to 
correspond to those of the infinite array analyzed by Amitay and Galindo (refs. 51 and 54). 
They employed a different time convention (e'^“^) in their analysis; therefore, a change 
in sign for their infinite- array reflection coefficient phase calculations (N =«>) was 
necessary for a direct comparison with the finite- array results. 

Data are presented for the two finite- array sizes indicated by the dashed circles 
inscribed on the array grid in figure 17. The elements in each finite array (N = 37 and 
N = 183) are those whose centers lie either on or inside the dashed circle. Data are 
presented for the center element (C) for both array sizes and for two edge elements 
(A and B) of the larger array. 

The data are presented as a fimction of the differential phase shifts and 
between elements in the H-plane and E-plane directions, respectively. These phase 
shifts are related to the beam-pointing directional cosines T^ and T^ by (ref. 54) 

= 2t7 (0,714) T^ (171) 

.Ay = 27 t (0.714) sin 60° Ty (172) 


The finite-phased- array calculations are obtained by first determining the complex 
scattering matrix for the array. Then the complex amplitudes of the incident waveguide 
fields (ap ) are given a phase differential according to 
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Figiire l6. Dimensions for equilateral triangular grid array 
of circular waveguide apertures excited in mode. 
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H-Plane 


Figure I7.- Finite-size phased array of circular waveguide 
apertures in an equilateral triangular grid arrangement. 
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(173) 





= e 


X 1 



In order to scan the beam in the E-plane, i/-„ is set equal to zero and is varied; 

X y 

and, in order to scan the beam in the H-plane, is set equal to zero and is varied. 

Simultaneous variation of and would scan the beam in some other direction as 

determined by the directional cosines T and T . At each value of </< and «/< , the 

ratio of b (determined from the product of the scattering matrix [S] and the column 
^i 

matrix [a]) to ap^ determines the amplitude and phase of the reflection coefficient of 

the ith aperture with all elements excited so as to point the beam in the direction specified 
by and T^. 

The amplitude and phase of the reflection coefficient for the center element of the 
two finite-size arrays is presented in figimes 18 and 19 together with the infinite- array 



0 30 60 90 120 150 180 


'“y , deg 

Figure l 8 .- TEn mode reflection coefficient of center element 
of arrays in figure I7 as a function of E -plane scan. 
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calculations (shown dashed) of Amitay and Galindo (refs. 51 and 54). As the size of the 
array is increased, the reflection coefficient exhibits a resonance behavior corresponding 
to the "blind spot" of the infinite array. Although the reflection coefficient of the finite 
array never reaches unity, the qualitative agreement with the infinite array calculations 
tends to justify the present analysis as applied to finite arrays. 



**'x . deg 


Figure I9.- mode reflection coefficient of center 

element of arrays in figure I7 as a function of 
H-plane scan. 


Figures 20 and 21 show a comparison between the reflection coefficients of the center 
element and the edge elements of the larger array (N = 183). Notice that the reflection 
coefficient of the edge element exhibits a much sharper resonance when the array is scanned 
in one direction. When the array is scanned in the opposite direction, the edge element 
reflection coefficient is almost constant and very near the isolated element value (N = 1), 
This asymmetry with scan is a general characteristic of the edge elements in a large 
periodic array. (See ref. 98.) 
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The resonance phenomena which occurs for the edge element when the array is 
scanned in one direction can be attributed to the destructive interference between the 
direct radiation from the edge element and a leaky wave traveling in one direction on the 
periodic structure (ref. 30). In the case of the center element, leaky waves traveling in 
both directions can produce symmetrical element pattern interference nulls (ref. 30) 
or impedance resonances. 



“'y , deg 


Figure 20.- TEn mode reflection coefficient as a function of 
E-plane scan for center and edge elements of l85-element 
array. 
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360 



Figure 21.- TE3_]_ mode reflection coefficient as a function of 
H-plane scan for center and edge elements of I83 -element 
array . 


The calculations in figures 22 to 24 are for the array geometry of figures 16 and 17 
with a dielectric cover of 0.5-wavelength thickness and dielectric constant of 2.0. In order 
to avoid numerical difficulties, a dielectric loss tangent of 0.0001 was assmned for the 
finite-array calculations. Previous results for the self-admittance of a dielectric-covered 
rectangular slot (ref. 124) indicated that a loss tangent of 0.001 or less would yield 3 or 4 
significant figure accuracy when compared with calculations for a lossless dielectric. 

In figure 22, the calculations for the center element reflection coefficient of the 
larger finite array (N = 183) exhibit two peaks which appear to correspond to the 
resonances of the infinite array (ref. 54). In order to verify this result, the reflection 
coefficient amplitude of the edge element B is compared in figure 23 with the infinite- 
array calculations. When the array is scanned in one direction, the edge element "sees" 
a much larger periodic structure and the destructive interference mentioned earlier pro- 
duces two sharper resonant peaks near the infinite -array "blind spots." This qualitative 
agreement between the dielectric- covered infinite- array and large finite- array calculations 
tends to establish further the validity of the present analysis. 
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Amplitude 



Figure 22.- TE^ mode reflection coefficient for the center 
element of arrays in figure I7 with a dielectric cover 
(dielectric constant, 2.0; loss tangent, 0.0001; dielec- 
tric thickness, O.5A; E-plane scan). 



Figure 25*- TEn mode reflection coefficient for edge element of 
183-element array with dielectric cover (dielectric constant, 2,0 
loss tangent, 0.0001; dielectric thickness, 0.5^i E-plane scan). 






'♦'x , deg 

Figure 2 k. - TE3_]_ mode reflection coefficient for the 
center element of arrays in figure I7 with a 
dielectric cover (dielectric constant, 2.0; loss 
tangent, 0.0001; dielectric thickness, 0.5'^; 
H-plane scan) . 


The reflection coefficient of the center element of the two dielectric- covered finite 
arrays is presented in figure 24 as a function of the H-plane scan parameter. The infinite- 
array calculations were not available for comparison. The finite- array calculations for 
this case are included since a different type of resonance occurs which probably warrants 
further investigation. It appears to be related to the phased- array impedance-matching 
techniques by the dielectric loading of the aperture plane (refs. 78 to 80) and indicates that 
the present analysis may also be of some benefit in the impedance matchir^ of finite 
arrays. Minor modifications (ref. 134) of the present analysis could also be used to study 
the impedance properties of finite arrays with dielectric plugs. 
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CONCLUDING REMARKS 


A variational expression has been derived for the self and mutual admittances of 
waveguide-fed apertures radiating into a multilayered region which may contain inhomo- 
geneous layers. The general expression has been evaluated for circular apertures 
excited in the TEmn TMmn waveguide modes and a computer program written 

which can include up to four external layers, two of which may be inhomogeneous normal 
to the aperture plane. 

Good agreement was obtained between measured and calculated values for the trans- 
verse electric (TEji) mode mutual coupling of two circular waveguides radiating into free 
space and one dielectric layer. A comparison was made between measured and calculated 
results for several combinations of frequency, polarization, and spacing. Very good 
agreement was obtained in all cases, except where the diffractions from the edges of the 
30.48 cm by 60.96 cm (12 in. by 24 in.) ground plane produced large scatter in the meas- 
ured data. 

By performing a parametric study, it was determined that the only significant effect 
of higher order modes is due to the transverse magnetic (TMn) mode and then primarily 
in the phase of the coupling coefficient. 

A comparison was also made between the reflection coefficient of an infinite array 
and the reflection coefficients of several elements of two finite arrays. It was shown that 
the center element of a 183-element array (approximately 10 wavelengths wide) had similar 
impedance characteristics to that of the infinite array. Although total reflection did not 
occur in the finite array, a definite peak in the reflection coefficient, corresponding to the 
"blind spot" of the infinite array, was observed. 

The validity of the theoretical model has been established by a comparison with 
measurements on two circular waveguide-fed apertures and with calculations on an 
infinite periodic array. 

This work allows the determination of mutual coupling between the elements of a 
finite array of circular apertures including any number of higher order modes and an 
arbitrary polarization for each aperture field. 

The analysis presented here can also be applied to other aperture shapes for which 
the Fourier transforms of the aperture electric fields can be determined. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., July 3, 1973. 
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APPENDIX A 


VECTOR POTENTIALS AND WAVE EQUATIONS 

The definition of the vector potentials (A and F) in their relationship to the 
electromagnetic fields (E and H) and the wave equations, as used in this analysis, 
are given here for reference. 

It is assumed throughout the paper that the electromagnetic fields contain a harmonic 
time variation of the form then. Maxwell's equations for a charge-free region are 


V X H - jcoe E = 0 

(Al) 

V X E + - 0 

(A2) 

V. yaii = 0 

(A3) 

V . e E = P 

(A4) 

where the permittivity e and permeability m may be complex and also be a function of 
the X., y., and coordinate variables. 

From equation (A3), the vector can be defined as the curl of another vector A, 

that is. 

ii =lv X A 

(A5) 

Substituting equation (A 5) into equation (A2) gives 


vx (E + jxuA) = 0 

(A6) 

Or since a curl-free vector is the gradient of a scalar 


E + j o)A = - V >P 

(A7) 

the electric field is given by 


E = - V»p - jcoA 

(A8) 


Similarly, by defining another vector 


F which satisfies equation (A4) so that 
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APPENDIX A - Continued 


E = - A V X F 

e 

(A9) 

then from equation (Al) 


H = - V$ - joiF 

(AlO) 

By superimposing the results due to the assumed magnetic vector and scalar potentials 
and the electric vector and scalar potentials, 

E = -AvxF-v>P-j CO A 

€ 

(All) 

H =iv X A - VO - jwF 

(A12) 


describes the electromagnetic fields in terms of a set of arbitrary vector and scalar 
functions. 

Since these functions are arbitrary, A is chosen as 


A = A Zj 


(A13) 


where A is a function of x., y., z.. Then with <i) and F temporarily set to zero 
and by using the vector identities 


V X 


vx A 


= V[l\x[VxA] +1[VXVXA] 


VXVXA^V(V. A) A 


(A14) 

(A15) 


equation (Al) becomes 


3 /1\ 3A ^ 1 3^A . 

+ 


3Z. BX. M BX-BZ. BXj 


Xi + 


3 /I 


BA ^ 1 B^A . BW 
+ + ]ue 


BZj \M/ At By^BZj^ 


B /1\ BA 


Bx. \/x/ Bx. 


B 


BA , IB^A . BW 2 . 1 * 

+ + Jcoe - oj eA-_V A 

Bz. fi 

'1 


^i 


^ Bz^ 


i. ^0 


(A16) 
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APPENDIX A - Continued 
Equating either the or y. components to zero gives 

'1'= - Jl (A17) 

Then the z. components yield 

And if the medium is assumed to be homogeneous in the and y. directions, equa- 
tion (A18) yields the wave equation 


1 Be B 
e BZj BZj 


+ k" 



(A18) 


v2 

A(xp y., zp 

_ , 1 3 

A(x.,y^z.) 

^^2 

A(Xj, y., zp 


1 

1 

e(z^ 3Z. BZj 

m(Zi) 

+ Kq 

^0^0 



(A19) 

Likewise, if it is assumed that "P = 0 and A = 0 momentarily and 

F = F z. (A20) 

then from equation (A2) 

<S>=-^±(1) (A21) 

iufj. BZ^ \€ / 

and, for /x =/x(z^) and e = e(zj). 


F(xj,yj, z.) 

1 B 

F(Xj,yj,Zjr 


mCz^) e (zp 

F(Xi, Yj, z^ 

e (Zi) 

/x(z.) Bz. Bz. 

e(Zi) J 

M e 

0 0 

e(Zi) 


(A22) 

Therefore, by superposition, the electric and magnetic fields can be derived from a 
set of z^ directed vector potentials which satisfy the differential equations (A19) and 
(A22) subject to the appropriate boimdary conditions on the electric and magnetic fields. 
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APPENDIX A - Continued 


If bidimensional Fourier transforms are assumed so that 


00 C 

A(k^,ky,Zj)= f j 

*Lr 


A(Xj.yj,z,) -jk^Xj -Jk yj 


00 •'-00 


e e ^ dxj dy. 


^ F(x.,y.,z.) -jk X. -jk y 


- r r - -jk X. -jk y. 

E(kjj,ky,Zi) = I j E(x.,y.,z.) e ^e ^ ^ dx- dy. 


•'-00 *'^co 


(SO jxi _ -jk„x. -jk y. 

H(k^, ky , ^i) = J J yi» Zj) e ^ e ^ ^ dx. dy. 

•^00 -00 


and inversely 


A ( xj , y ,, zj ) 1 


( 2»)2 


r“ r“ j^^i j^x^yi 

J J A(kjj.kj,,z^) e^ ‘e y *dk^ dky 


F(Xi.yi.Zi) I 


J / *x 

nn • L/vx 


( 2 t 7 )^ * 5-00 •'-00 


dk. 


jk X. jk y. 
J X 1 J y-^l 


'i»yp^i)=-^ f f ^ ^dk^dk^ 

(277^ •’-00 •'-00 


1 r" r“ - j^v^i j^x^yi 

H(Xj.yj,z,)= j f H(k^,ky,Zj)e eF'dk^dkj 

( 2 t 7 ) *^00 •'-00 


(A23) 


(A24) 


(A25) 


(A26) 


(A27) 


(A28) 


(A29) 


(A30) 
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APPENDIX A - Continued 


then the transverse components of the transformed electric and magnetic fields become 


1 


(A31) 




yi ^ y ^ coeiz.) 


(A32) 


Hx,0Sc> S’ ‘ V 


(A33) 


Hy,OSc.Vl>'l‘ScA»Sc’V^>-;^) 


f'0sc> S’ ^i> 


(A34) 


where the primes on A and F denote differentiation with respect to z^, and where 
the transformed potential functions now satisfy the differential equations 


^2 1 de(z.) ^ 

^ ^ AA(k^, ky, z.) 


A(kx» k , z.) - 

2 ^ y e(zj) dzj dZj^ 


1 2 
+ kp 


M' 


(.i) 4^1) 4 + 


f^n^i 


0 0 


A(k^, k , zp = 0 


(A3 5) 


« F(k^,k k;)- ' 


dM(Zi) 


dz^ 


2 ^ y i dz. dz. 


Fds,, k z.) 




IJ.{Z^) e(zp + ky 


^0^0 


"■0 


F(kjj, k , Zj) = 0 


(A36) 
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APPENDIX A - Concluded 


At a boundary (z^ = dp) between two media (a^(zP €p(Zj) and Mp+i(Zi) ep+i(Zi)), 
continuity of the transverse electric and magnetic fields or their transforms shows that 
the transformed potentials for the two regions must satisfy the boundary conditions 


ApOSt. ky, dp) = Ap^j(k^, ky. dp) 

PpAsc’ is,’ V = ‘V- V 


(A37) 

(A38) 


^VscV“i> 


)^V r 








dz. 


z.=d 
1 P 


(A39) 


Tif S’ 


. dz 
L 1 


-Izpdp 


57 N’ 

L^i 


^i=dp 


(A40) 
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APPENDIX B 


COMPUTER PROGRAM FOR THE CALCULATION OF THE SCATTERING MATRIX 
OF A PLANAR ARRAY OF CIRCULAR WAVEGUIDES RADIATING 
INTO EITHER FREE SPACE OR FOUR DIELECTRIC LAYERS 

The program listed here computes the complex mutual admittances between the modal 
fields of all the apertures in the array. These admittance values are used to form a com- 
plex square matrix which is operated on with the appropriate matrix algebra and inversion 
to obtain the complex scattering matrix for the array. 

The basic program for the mutual admittance calculation is a modification of a 
previous computer program (ref. 135) for the calculation of the self admittance of a TEjj 
mode excited circular aperture. The computer program in its present state is limited 
to a maximum of four external layers over the apertures; however, the third and fourth 
layers may be inhomogeneous normal to the apertime plane. 

The wave equations (eqs. (58) and (59) with m(z.) = Constant) for the third and fourth 
layers are solved numerically by using a spline routine for curve fitting through discrete 
points of the dielectric or plasma profile. The numerical integration of equations (167) to 
(170) is performed by a Runge-Kutta method containing a variable increment which is 
continuously subdivided until a specified :accuracy is achieved over each integration step. 
More detailed discussion of the numerical techniques used in the original computer pro- 
gram are given in reference 135. 

The present program accepts as input the following parameters; 

NUMMODE total number of waveguide modes assumed in each aperture (each 

aperture field distribution is assumed to be a superposition of the 
same waveguide modes) 

NUMHOLE total number of circular apertures in the array 

NUMTE total number of transverse electric waveguide modes assumed in 

each aperture field distribution 

NUMTM total number of transverse magnetic waveguide modes assiuned in 

each aperture field distribution 

(Note that NUMMODE = NUMTE + NUMTM) 

MIJ(I), NIJ(I) indices of TE^ modes, 1 = 1 to NUMTE (if NUMTE = 0, omit) 
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APPENDIX B - Continued 

MMIJ(I), NNIJ(I) indices of modes, I = 1 to NUMTM (if NUMTM = 0, omit) 

AIJ(I) radius (a^ of each aperture, I = 1 to NUMHOLE 

XI(I), YI(I) x,y coordinates of center of each aperture (x!, y|), I = 1 to NUMHOLE 

PHIJP(I) angular rotation (</>p of x^ axis with respect to x axis for each 

aperture, 1=1 to NUMHOLE 

F frequency 

ZI, Z2, Z3, Z4 distances (dj, ^ 2 , dg, d^) from aperture plane to outer surfaces 

of layers 1, 2, 3, 4, respectively (For radiation into free space, 
set Z4 = 0.0.) 

CONVERT a conversion factor to change all input dimensions to centimeters 

(that is, for input dimensions in inches, set CONVERT = 2.54) 

ER relative dielectric constant of material completely filling all wave- 

guides (ER = 1.0 for air-filled guides) 

For radiation into free space (Z4 = 0.0), no additional input data is needed; however, 
for Z4 > 0.0, the relative values of Zl, Z2, Z3, and Z4 are compared to determine 
which external layers are to be considered in the calculations and appropriate parameters 
are read in as follows: 

Vl,Wl real and imaginary parts of complex relative dielectric constant of 

layer nearest to aperture plane (layer 1) (if Z3 = 0.0, and Z4 > 0.0, 
omit) 

V2,W2 real and imaginary parts of complex relative dielectric constant of 

layer 2 (if Z3 = 0.0 and Z4 > 0.0, omit) 

NP3 number of points used in approximation of the inhomogeneous profile 

for the dielectric constant of layer 3 (if Z3 = 0.0 and Z4 > 0.0, 
omit; or if Z2 = Z3 = Z4, omit) 

ZD(I) distance from aperture plane to discrete points in layer 3 dielectric 

profile (if Z3 = 0.0 and Z4 > 0.0, or if Z2 = Z3 = Z4, omit) 

V3(I), W3(I) real and imaginary parts of relative dielectric constant at discrete 

points ZD(I) in layer 3 inhomogeneous profile (if Z3 = 0.0 and 
Z4 > 0.0, or if Z2 = Z3 = Z4, omit) 
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APPENDIX B - Continued 


NP4 number of points used in approximation of the inhomogeneous profile 

for the electron density and collision frequency for the plasma of 
layer 4 (if Z3 = Z4, omit) 

ZND(I) distance from aperture plane to discrete points in layer 4 plasma 

profile (if Z3 = Z4, omit) 


NE(I), NU(I) electron density and collision frequency at discrete points, ZND(I), in 

layer 3 inhomogeneous plasma profile (if Z3 = Z4, omit) 

The values of ZD(I) and ZND(I) must be monotonically increasing. The first 
value of ZD(I) must equal Z2, the last value of ZD (I) and the first value of ZND(I) 
must equal Z3, and the last value of ZND(I) must equal Z4. Any deviation from these 
values will cause errors to occvm in the calculations. 

The output of the computer program is as follows: 

YC(II,JJ) 

PRMT(2) 

IHOLE 
JHOLE 
IMODE 
JMODE 
YMN(I) 

S(I, J) 


elements of complex admittance matrix 

' ! 

upper limit of numerical integration (maximum value is set at 50.0) 
ith aperture , 

jth aperture 

pth mode in ith aperture 
qth mode in jth aperture 

characteristic admittance of pth mode in ith aperture 
elements of complex scattering matrix 
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APPENDIX B - Continued 


PROGRAM CIRWG( I NPUT • OUTPUT • TAPES* I NPUT • TAPE6-0UTPUT ) 


C 

C 

C 

C 

c 

c 


MI •Nl ARE 
AT = RADIUS OF 
RADIUS OF 
CENTER TO 


AJ 

R 

PHI 
PH I JPP 


the subscripts 
I-TH aperture 
U-TH APERTURE 

center spacing 


OF THE aperture 


h *»««»«» «**««*««« 

excitation modes 


S ANGULAR rotation OF R WITH RESPECT TO XI-AXIS (DEGREES) 

= ANGULAR rotation OF XJ-AXIS WITH RESPECT TO XI-AXIS (DEG.) 

•*******♦•»*»*#**♦ 


D I MENS ION AUX ( 8 « 4 ) . DpRY ( 4 ) « PRMT ( 5 ) » Y ( 4 ) 

O I MENS ION ZD ( 50 ) ♦ ZNO ( 50 ) ♦ NE ( 50 ) « NU ( 50 ) 

C DIMENSIONS FOR COMMON VARIABLES 

DIMENSION V3 (50 ) . V31 (50 ) . V32 ( 50 ) . V33 ( 50 ) . 

A V4 (50 ) « V41 (50 ) ♦ V42 ( 50 ) ♦ V43 ( 50 ) ♦ 

B W4 (50 ) »W41 (50 ) . W42 (50 ) »W4 3 (50 ) * 

C W3(50 ) .W31 (50 ) .W32 (50 ) «W33 (50 ) . 

D Z(50 ) .ZN(5n ) ,XMNP(8.3 ) «XMN (8* 3 ) 

C COMMON - DIMENSIONED VARIABLES 

COMMON V3.V31 ♦V32*V33*W3«W31 «W32 »W33 »Z«XMNP *XMN* 

A V4.V41 »V42*V43»W4*W41 «W42»W43«ZN* 

C COMMON - UNoIMENSI ONEO VARIABLES 

B BSO. CCA. CCB.D?. KINO. L3.L4, MOST. TMI .TMJ, 

C NEW. NP3. NP4. RkERR. RK3. RK4 .TERMA, TERMS, TERMC.TERMD , 

0 VI. V2. V3XI3.V4XI4. 

E W1 .W1SQ.W2.W2S0.W3XI3.W4XI4.XI 1.XI2.XI3.X14, 

F Ml .MJ.Nl .NJ.MIP.MJP,AKZER01 . AKZEROJ . RKZERO . FACTOR I .FACTORJ. 

G FACtSQI .FACTSOJ.COSP.COSM.SINP. SINM .COSPHIJ.PHI JPP 
LOGICAL TMI.TMJ 

REAL KZERO.NE.NU.NED.NUD. ISOLATE 

complex CCA . CCB . COEFF . GAMMA . CON . YC . YTE . YTM , YMNZERO 
complex a.b. determ. yap 

01 MENS 1 ON YC (25.25) . A { 25 .25 ) .B (25 .25 ) . I p I VOT ( 25 ) . I NOEX (25.2) 
dimension A1J(25).XI (25).YI (25).PHIJP{25).MIJC10),NIJ(10) 

0 1 MENS I ON MM IJ ( 1 0 ) . NN I J ( 1 0 ) 

EXTERNAL FINDC 

C fSTABLISH CONSTANTS 

XMNP (1.1 ) =3.832S>XMNP (1.2) =7. 01 6SXMNP ( 1 , 3 ) = 1 0. 1 73SXMNP (2.1 ) = 1.84118 
XMNP (2. 2 ) ®5.331 tXMNP (2.3 ) =8.536SXMNP (3.1) =3.054SXMNP (3.2)=6.706 
XMNP (3.3) =9»969iXMNP (4.1 ) =4 .201 SxMNP {4,2;=8»015SXMNP(5.1 )=5.317 
XMNP (5. 2 ) =9.282*XMNP (6 , 1 ) =6.4 1 6SXMNP (6 . 2 ) = 1 0.52SXMNP (7,1 ) =7. 501 
XMNP (8. 1 )=8.578 

XMN(1 .1 )=2.405SXMN(2 , 1 ) = 3,832$XMN(3, 1 ) =5 . 1 36T.X MN ( 1 ,2)=5.520 
XMN(4.1 )=6.380*XMN(2,2)=7,016$XMN(5. I ) =7 . 588SX MN ( 3 . 2 ) =8 . 4 1 7 
XMN( 1 . 3 )=8.654SXMN(6 , 1 ) =8 . 77 I $XMN ( 4 , 2 ) =9 , 76 1 $XMN ( 7 , 1 ) =9,9 36 
XMN(2.3 ) = 10, 1 73SXMN (F. 2 ) = 1) .065SXMN ( 8. 1 ) = 1 1 . 086 
P1=2.0*AS1N( 1 .0 ) 

TW0PI=2.0*PI 

F0RK = TW0PI/(3.*1 .El O ) 

FOMEGAs (TWOPI*8970, )#*2 
C0N=cMPLX(0.0. 1 .0 ) 

C START READING INduT 

C NUMMODE = NUMBER OF MODES PER APERTURE 

C NljMHOLE = NUMBER OF APERTURES 

C NUMTE = NUMBER OF TE MODES 

C NuMTM = NUMBER OF TM MODES 

1 read (5. 1 1 )NUMH0LE. NUMMODE, NUMTE. NUMTM 
11 format (4 15) 
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lF<e0F«5)900«2 

2 M=NUMHOLE*NUMMODe 
tF(M,GT.25)3.4 

3 WRITf(6»80) 

80 format (1 HI *M EXCEEDS DIMENSION OF YC* ) 

STOP 

4 IF(NUMTE«GT.lO)GO TO 5 
IFINuMTM.GT. 10 )G0 TO 5 
IFINUMTE.GT.NUMMODEIGO to 900 
IF (NUMTM.GT.NUMMODE )G0 TO 900 
IF <NUMMOOE.GT.20 )5.6 

5 WRITf(6«81 ) 

81 format ( 1 HI *NUMM0DE exceeds DIMENSION OF MU AND N IJ* ) 

STOP 

6 IF <NUMH0LE«GT«25 )7. 1 F 

7 WR1TF(6«84) 

84 format ( 1 HI »NUMH0LE EXCEEDS DIMENSION OF A IJ* ) 

STOP 

C*********************************************************************** 
c MIJ(I),N1J(I) = INDICES OF I -TH MODE TE-MN 

C MMlJ( 1 )«NNIJ< I ) = INDICES OF I-TH TM-MN MODE 

c******************************** ****************■"■*■*■ *********** ********* 
12 IF(NuMTE.EQ.O)GO TO 17 

read (5«14)((MIJ(I)«NIJ(I))«1=1« NUMTE ) 

17 IF (NUMTM.EQ.O )GO TO lO 

REAO(5« 14) ( (MM1J( I ) «NNIJ( I ) ) « 1 = 1- .NUMTM ) 

10 CONTINUE 

14 F0RMaT(20(2I1.2X)) 

IF (NUMTE. EO.O )G0 TO 
DO 50 1=1 .NUMTE 

IF (Ml J ( I ) .GT.7.0R.N I j ( I ) »GT. 3 ) GO TO 601 

50 CONTINUE 

55 IF (NUMTM.EQ.O )GO TO 53 

DO 52 1=1. NUMTM ' 

IF (MM I J ( I ) .GT.7.0R.NNI J( I ) .GT.3 )G0 TO 602 \ 

52 CONTINUE 

53 CONTINUE 

C*********** *************************************************** ********* 

c AIJ(I) = radius of i-th aperture 

C XI (I) AND YI(I) = X.Y COORDINATES OF CENTER OF I-TH APERTURE 

C PHIJP(l) = ANGULAR ROTATION OF XI -AX IS WITH RESPECT TO X-AXlS 

C (DEGREES COUNTER-CLOCKWISE). 

^************ **»**»**♦*#♦***♦**♦ #»*************»**♦♦***** ****** ********* 
READ(5.15) ( A I J ( I ) . 1 =I .NUMHOLE ) 

READ(5.15) (XI ( I ) .1 = 1 .NUMHOLE) 
read (5. 15) (Y I ( I ) . 1 =1 .NUMHOLE ) 

RE AO ( 5 . 1 5 ) ( PH IJP ( 1 ) . I = 1 . NUMHOLE ) 

15 FORM aT(8F1 0.2 ) 

S1ZE=1 .0 
P0L=0.0 

IF (NUMHOLE. GT. 1 )8«9 

8 DO 60 I =2. NUMHOLE 
J=I-1 

IF(AIJ( I ).NE.AI J( J) ) S1ZE=0.0 

60 CONTINUE 

DO 51 I =2. NUMHOLE 
J= 1-t 

IF(PHl JP( I ) .NE.PHIJP( J) ) P0L=1.0 

51 CONTINUE 

9 CONTINUE 
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C F = frequency ( CYCLES PER SECOND) 

C Z1«Z?*Z3«Z4 = distances FROM APERTURE TO OUTER SURFACE OF LAYERS 

C 1 ♦2.3.4 RESPECTIVELY. 

C convert = CONVERSION FACTOR FOR CONVERTING INPUT DIMENSION TO 

C centimeters (IF INPUT IN CENTIMETERS. LEAVE BLANK), 

C ER = relative dielectric CONSTANT OF MATERIAL FILLING ALL WAVEGUIDES 

(;»*♦*#♦*#*♦*******»**♦*#***♦***♦*■»*«*«•***■*■****♦*♦♦***■»**»**»#■*#*♦*♦♦#*»* 

C FOR all input dimensions in wavelengths 

C set F*3.0E10 and CONVERT=1.0 

c 

c FOR FREE SPACE. SET Z4=0.0 

C 

READ(5. 16)F,Zt . 22 . Z3 , Z4 , CONVERT . ER 
IF(ER.EO.O,0 )FR=1 ,n 
16 F0RMAT(E10,2.7F10,2) 

WRITc-(6.41 ) 

41 FORMAT! 1 HI) 

RKERR=0,OOOl 

EPSLN=0.00001 

C RKERR.EPSLN = ERROR TOLERANCES FOR RUNGE-KUTTA INTEGRATION 

C AND SPLINE CURVE-FIT ROUTINE. 

DIMEN=1 .0 

IF(CONVERT.NE,0.0 )DIMEN=C0NVERT 
WRITF(6.105) 

IF(Z4.NE,0.0 )G0 TO lo 
WR I TF ( 6 . 1 8 ) F . D I MEN 

18 format ( 3X*MUTUAL COUPLING OF CIRCULAR APERTURES RADIATING INTO FRE 
IE SPACE*//5X#F = *E1?.5/5X»DIMEN = ♦F10.6/) 

KIND=4 
VI =V2=1 .0 
W1“W2*0»0 
GO TO 21 

19 CONTINUE 

WRITf<6.20 )F.Z1 .Z2.Z3.Z4.0IMEN 

20 format (3X4MUTUAL COUPLING OF CIRCULAR APERTURES RADIATING INTO A M 
IULTIlAYERED dielectric under a NONHOMOGENEOUS plasma LAyER.*// 5X» 
2INPUTS»//2X»F a *E12,5/1X*Z1 = *F10.6/1X*Z2 = »F10.6/1X*Z3 = *F10. 
36/1X4Z4 = *F10,6//1X*DIMEN = *F10.6) 

WRITf(6. 1 05) 

21 CONTINUE 
WR1Tf(6.85)ER 

85 format (1X»ER = *F10.6) 

DO 70 1=1.NUMH0LE 

WRITe< 6.92 ) 1 .AI J ( 1 ) . I .XI ( 1 ) . I .Yl ( 1 ) . I .PHI JP( 1 ) 

92 FORMAT! 1 X*A I J ( # 1 2# ) =#F8 . 5 . 5X*X I ! * I 2* ) =*F8. 5 . 3X*Y I ! * 1 2* ) =*F8. 5 .5X*P 
1HIJP(*I2*)=»F8.3* DEG.*) 

70 CONTINUE 

93 FORMAT! 1X*M00E*I2* = TE-*2I1) 

WRITr(6. 1 05) 

DO 72 I=1.NUMM0DE 

IF! I .GT,NUMTE)GO TO 7l 

WR I TF ( 6 . 93 ) I . M I J ( I ) ♦ N I J ! I ) 

GO TO 72 

71 IOEM=I-NUMTE 

WRITf!6 .91 ) I .MMI J ! I OfM ) ,NNI J ! IDEM ) 

91 FORMAT! 1X*M0DE» I 2» = TM-*2I1) 

72 CONTINUE 
WRITf (6.105) 

IF!Z4.EQ.0.0 )G0 TO inO 
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C PIND INPUT DATA rASP 

KIND=1 

1F(Z4-23)26*47,33 

26 WRITf(6«27) 

27 FORMAT! 1 OX* 16H INPUT DATA ERROR) 
GO TO 900 

33 IF (Z3 )26*40«82 
40 KIND=3 



GO TO 

1 09 

47 

IF !Z7- 

Z2 I26f54«61 

54 

KINDsA 



GO TO 

82 

61 

K1NDs2 



C GET OTHER INPUT NEEDED BASED ON INPUT DATA CASE 

C VI. W1 = REAL AND IMAGINARY PARTS OF DIELECTRIC CONSTANT OF LAYER 1. 

C V2«W2 = real and IMAGINARY PARTS OF DIELECTRIC CONSTANT OF LAYER 2. 

82 REA0(5. 13)V1 *W1 «V2.W? 

WRITP(6.105) 

WRITf(6*83)V1 *W1 *V2.W2 

83 F0RMAT(3X«V1 = *F8.5,5X*W1 = *F8.5/3X#V2 = *F8.5.5X»W2 = »F8.5/) 

. IF(KIND.EQ.4 ) GO TO 1 00 

««««*««««*»****************«'«'**«»****«««"»**)!'*!)))!*)!* *-******)))!** )))^ *«)!))*** 
C NP3 = NUMBER OF POINTS FOR LAYER 3 DIELECTRIC PROFILE 

(.*********** »*«**#*.|HHHf*##*»**»»**********#*»*»*#)H(-************#)HHHHMH»-* 

READ(5* 1 I0)NP3 
IF (NP3 )89»89*96 

89 WRITF(6«90) 

90 FORMAT! 1 Ox. 42HERR0R IN NUMBER OF POINTS FOR V3 W3 TABLES/) 

GO TO 900 

96 IF !NP3-50 ) 1 02. 1 02 .89 ' 

^****#******»*##*#»**»**##******#*«*****##*****»#****#*4»»*-)HHH!#*#*-)H!#*# 
C V3!I).W!3) = REAL AND IMAGINARY PARTS OF DIELECTRIC CONSTANT AT 
C POINTS Zn!I) INSIDE LAYER 3. 

102 READ !5. 15) ! ZD ! 1 ) . I = 1 . NP3 ) 

READ!5.15) !V3! I ) . I = 1 .NP3) 

READ! 5. 15) ! W3 ! ! ) . 1 = 1 .NP3 ) 

WRITF!6.103) 

103 format ! 14X. 1 HZ. 15X.5hV3!Z ) . 1 2X . 5HW3 ! Z ) / ) 

WRITE (6. 1 04 ) !ZD ! I ) . V3 ! I ) .W3 ! 1 ) . I = I .NP3 ) 

104 F0RMAT!5X.3F17.5) 
write (6. 105) 

1 05 format ! 1H ) 

IF!K1N0»EQ*2) GO TO 1 00 

C '••"X' * X- ***)*"*'*'*"))"*"*'**'*')*"*"****'*******)! * 

C NP4 s NUMBER OF POINTS FOR LAYER 4 PLASMA PROFILE 

109 READ!5. 1 1 0)NP4 

1 10 FORMAT! 15) 

IF !NP4 ) 1 16* I 16* 123 

116 WRITf!6.1 17) 

117 FORMAT! 1 OX. 42HERR0R IN NUMBER OF POINTS FOR NE NU TABLES/) 

GO TO 900 

123 IF !NP4-50 ) 1 30. 1 30. 1 1 6 

(.**♦*****♦*♦»♦♦♦»*♦#♦*»##*»*»*****#*****#**•»***■*#♦#**♦♦♦*********♦»***♦* 
C NE!I).NU!1) s electron DENSITY !ELECTRONS PER CUBIC CENT I METER ) AND 
C electron collision frequency !PER SECOND) AT POINTS 

C ZND!I) INSIDE LAYER 4. 
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c**#********»***************-************ 

130 READ(5«15) ( ZNO ( I ) « I = 1 »NP4 ) 

REA0(5«13) (NE( t ) ♦ I *1 «NP4 ) 

REA0{5.13) (NU( I ) « 1*1 «NP4 ) 

13 format (BE 10, 2 ) 

WRITF(6«138) 

138 FORMAT! 14X» 1 HZ, 12X,5hNE (Z ) , 9X , 5HNU f Z ) / ) 

WR I Te < 6 , 1 39 ) ( ZNO ( 1 ) , NE ( I ) , NU ( 1) , I = 1 , NP4 ) 

139 F0RMAT(8X, 3E14.5) 

WRITF(6,105) 

100 CONTINUE 

WR1Tf(6,94 ) 

DO 1000 1H0LE=1 ,NUMHOLE 

DO lOOO JH0LE=1 ,NUMHOLE 

00 lOOO 1M0DE»1 »NUMM00E 
DO lOOO JM00E*1 ,NUMMODE 
TM1», FALSE, 

TMJ*. FALSE, 

1 F ( I MODE ,GT, NUMTE ) TM I » , TRUE , 

1 F ( JmODE , GT , NUMTE ) TM j= , TRUE , 

PHl=n,0 

1 1 * ( I HOLE- 1 ) »NUMMOOE+ 1 MODE 
JJ*<JH0LE-1 )*NUMMOOE+ JMODE 
IF ( 1h 0LE,E0, JHOLE )50o,503 

500 1F(SIZE,E0,1 ,0)501 ,503 

501 IF ( IHOLE.GT, 1 )502 ,503 

502 YC( I I , JJ)*YC( IMODE, JMODE) 

GO TO lOOO 

503 CONTINUE 

94 FORMAT! 1 HO* */ » 

A1=AI J! I HOLE )*D I men 
A J* A I J ! JHOLE ) *0 1 MEN 
IF!TmI 1504,505 

504 I OEMs IMODE -NUMTE 
PHPIsPHlJP! 1H0LE)*PI/180, 

Ml =MMI J ! IDEM ) 

N1 =NNl J ! IDEM ) 

GO TO 506 

505 MI=MIJ! IMODE) 

NI=N| J! IMODE) 

PHPIsPHlJP! IH0LE)*P1/180. 

506 IF!TmJ)507,508 

507 IDEMsJMODE-NUMTE 

PHP JsPH I JP ! JHOLE ) *P I / 1 80 , 

MJ=MMI J! IDEM) 

NJ=NIMI J( IDEM ) 

GO TO 509 ' 

508 MJ=MIJ! JMODE) 

NJ*NIJ! JMODE) 

PHPJsPHI JP! JHOLE )*P 1/1 80, \ 

509 CONTINUE 

IF!TMl ,AND, !MJ,EO,0 ) ,AND, I ,NOT,TMJ ) )G0 TO 99999 
IF!TmJ,AND, !MI ,EQ,0 ) ,AND, ! ,NOT,TMI ) )G0 TO 99999 
PH 1 JDP= PHP J-PHP I 
MIP=Ml+l ' 

MJP=MJ+1 

XJI=DIMEN*!XI ! JHOLE ) -XI ! I HOLE ) ) 

YJI=DIMEN*!YI ! JHOLE) -Y1 ! IHOLE ) ) 

R=SOPT!XJl*XJI+YJI*YjI ) 

IF ! IhOLE,EO, JHOLF )R = n,0 
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1F(AbS(PHIJPP).LT. 1 .OE-04 )PHI jpo=0.0 

lF{T«dI ♦AND# •NOT»TMJ)PHIJPP=PHIJPP-P I /2. 

IF<TMJ*AN0««N0T#TMI )PHIJPP=PH1 JPP-PI/2. 

1F(R.EQ»0.0). GO TO 703 
IF(AbS(XJ1 )»LT,1«0E-06)700*701 

700 PHI=0«5*P1 

IF(YJI .LT.O.O )PHt =PH1+Pl 

PHI=PHI-PHP1 

GO TO 70a. 

701 PHI=ATAN2(YJI «XJ!)-PhPI 

702 ARG* (MJ+MI )*PH|-MJ*PHIJPP 
COSP=COS(ARG) 

SINPsSlN(ARG) 

ARGo ( M J-M I ) *PH I -MJ*Ph I JPP 
COSMsCOS (ARG) 

SINM=SI N (ARG ) 

703 COSPh1J“COS(PHIJPP) 

kzero*fork*f 

A<ZEROI=AI*KZERO 

AKZER0J=AJ*KZER0 

rkzepo=r*kzero 

IF(Tm1 )704«705 

704 FACToRI=XMN(M!P.NI)/AKZEROI 
GO TO 706 

705 FACToRI =XMNP (MIP»NI ) /AKZEROI 

706 IF(TMJ)707*Tp6 

707 FACToRJ=XMN(MJP«NJ )/AKZER0J 
GO TO 709 

708 FACToRJ*XMNP ( MJP«NJ ) /AKZEROJ 

709 CONTINUE 
FACTb01=FACTORI**2 
FACTbOJsFACTORJ**2 
XI 1=71 »KZERd*OIMEN 
X I 2=72*KZFR0*D I MFN 
X I 3*73»KZER0*DI MEN 
X 1 4 = 74 »KZERO*D I MEN 
02=XT2-XI 1 
RK3=RKERR 
PK4=RKERR 

IF (R.EQ*0«0» AND»M 1 .NF»MJ ) GO TO 99999 

IF (KIND»EQ.3 ) GO TO 220 

W2S0=W2»W2 

CCA = CMPLX (VI «-Wl ) /CMPLX ( V2 *-W2 ) 

IF <K IND»NE«4 ) GO TO 210 
TERMC=V2 
TERM0=-W2 
GO TO 149 
210 CONTINUE v, 

DO 107 1»1»NP3 
Z( I )=ZO( I )*DIMEN 

107 2(11=2(1 14KZERO 

: 4Et up arrays for spline interpolation 

CALL SPLRE0(NP3«EPSLN«2»V3«V31 «V32*V33) 
call SPLRED (NP3»EPSLN«Z ♦ W3« W31 ♦ W32* W33 1 
L3=NP3 

CALL SPLD2 ( NP3 »L3»X I3«Z»V3«W3»V31 «W31 *V32»W32tV33«W33«V3X13» 
A W3X I 3« dummy. DUMMY ) 

IF(K1ND»EQ.2)G0 to 146 
220 CONTINUE 

pSTABLISH V4 and W4 


85 



APPENDIX B - Continued 


omega*twopi*f 

OMeG<^Q=OMEGA*OMEGA 
00 1*^7 Is1,NP4 
ZN( J )=ZN0( t )*0!MEN 
OPSOsFOMEGA»NE ( I ) 

0EN0m=0MEGSQ+NU( I )**? 

V4 ( I ) = 1 .-OPSQ/DENOM 

W4 ( I 1 =NU ( I ) *OPSO/ (OMfG A*DENOM ) 

137 CONTINUE 

DO 140 I = 1 «NP4 

140 ZN( 1 >=ZN( n*KZEPO 

CALL SPLRE0(NP4*EPSLN«ZN.V4 .V41 «V4a* V43 ) 

CALL SPLRE0(NP4«EPSLN«ZN.W4,W41 «W42«W43) 

U4SNP4 

CALL SPLD2(NP4*L4»XU*ZN«V4«W4«V41 .W41 . V4 2 . k'42 ♦ V4 3 . W43« VAX 1 4 « 

A W4X14,0UMMY,DUMMY 1 

CALL SPLD2 (NP4 .L4«XI3«ZN«V4 .W4. V4 1 , W41 . V42 ♦ W42 . V43 » W43 . V4 X 1 3 « 

A W4XI3»DUMmY.DUMMY ) 

<^ET UP INITIAL CONDITIONS FOR BASE RUNGE-KUTTA INTEGRATION 
IF(KlND»NE.l )G0 TO 148 
0EN0 m=V4X I 3**2+W4XI 3*#2 

TERMA= ( V3X I 3*V4X 1 3+W3X I 3*W4X I 3 ) /DENOM 
TERMRs ( V 3X I 3»W4X 1 3-W3X 1 3* VAX I 3 ) /DENOM 
146 L3=l 

CALL SPLD2(NP3«L3.XI2*Z» V3»W3.V31 ,W31 ♦ V32 « w32 . V33 . W33 , V3X 1 2 , 

A W3X12«DUMM Yt DUMMY ) 

DEN0M = V3X 1 2*»2+W3X I 2**2 
TERMC= ( V2*V3X I 2+W2* W3X 1 2 1 /DENOM 
TERMD* < V2*W3X 1 2-W24V3X I 2 ) /DENOM 
GO TO 149 

148 V1*V4XI3 
W1 =W4XI3 

149 W1SQ=W1*W1 
CCB=cMPLX(-V1 «W1 1 
PPMT( 1 )=0. 

NPRm = 0 

PRMT(2)=0.0l 

PRMT (3 1 = (PRMT (2 1-PRMT ( 1 ) )/5. 

Y1 TEST=Y2TEST=Y3TEST=Y4TEST=0.0 
PRMT (4 ) =RKERR 

150 PRMT (5 1=0, 

NEW = n 
MOSTcO 

DO 151 I=l«4 

Y(M=0, 

0ERY(I)=>.25 

151 CONTINUE 

CALL LRKSl ( PRMT ♦ Y «OERY , 4 « F I NDC* AUX ) 

IF (PRMT (5) It «165«158 
158 IF (NPRM»GE.2 )G0 TO I6l 

IF(PRMT(1 ).NE.O.O)GO TO 162 
WR 1 TE ( 6 « 1 57 ) PRMT ( 1 ) * PRMT ( 2 ) 

PRMT ( 1 ) =0,001 

PRMT ( 3 ) = ( PRMT ( 2 ) -ORMt ( 1 1 ) /5 • 

PRMT (4 ) =RKERR 
NPRM=0 
GO TO 150 

162 write (6« 1 57 )PRMT ( 1) ,ORMT(2 ) 

NPRM=NPRM+1 

PRMT ( 3 ) «PRMT ( 3 ) /5 » 



APPENDIX B - Continued 


GO TA 150 

16J prmT (4 ) =PRMT ( 4 ) * 1 0, 

WR I TF ( 6 ♦ 1 59 ) PRMT ( 4 ) 

159 FORMAT! 1 Ox. 42HERROR TOLERANCE 

A 10X.18HBEEN INCREASED TO .E12.5) 

WR I TF ( 6 . 1 57 ) PRMT { 1 ) . pRMT ( 2 ) 

J57 format ( 1 0X*PRMT ( 1 )=*F6.3.3X*PRMT (2 )=*E6.3 ) 

GO TO ISO 


FOR LRKSI INTEGRATION HAS ./ 


integration complete, calculate final answers 
Y d) = real part of yte integral 

V(2) = REAL PART OF YTM INTEGRAL 
Y(3) * IMAGINARY PART OF YTE INTEGRAL 

Y<4) = IMAGINARY PART OF YTM INTEGRAL 

165 CONTINUE 

I F (PRMT (2)*LT«6*G)1 66 *171 

166 Y1TE!^T = Y1TE5T+Y(1 > 

Y2TE^T= Y2TESt + Y (2 ) 

Y3TE«^T = Y3TE5T + Y ( 3 ) 


Y4TE^T = YATES T + Y (4 ) 
PRMT ( 1 ) =PRMT ( 2 ) 


NPRM=0 

IF ( ArS (PRMT (2)-1*0).lE» 1 tE-OSlRRMT (2 ) = 1 


• 0 


PRMT ( 1 ) =PRMT ( 2 ) 

if(abs(prmt( n-1 •0 )*le.i 


E-05)PRMT( 1 )=1 .00001 


PRDEL“0 .24 

IP (PRMT (2 ) .GE.O .25. AnO.PRMT (2 ) •LT.2.0 )PRDEL=0.2S 

1P(PrMT (2). GE. 2.0. AND. PRMT! 21.lt ,4 . O ) PRDEL=0 . 5 

IP (PRMT (2 1 .GE.4.0 )PRnEL=l .0 
PRMT|2)=PRMT!2)+PRDFL 

IF1AbS!PRMT!2)-1.0>.lE.1.E-05)PRMT!2)=0.99999 

PRMT ! 3 > = ( PRMT ( 2 ) -PRMT ! I 1 1 /5 . 


PRMT ! 4 1 =RKERR 
GO TO ISO 

171 IFIPRMT (2 ).GE.50,0)G0 TO 175 
IFIAbSIYITEST) .LT. l .F-2001G0 TO 172 
1 F!AbS!Y( 1 )/Yl TEST 1-1 .E-04) 172.166. 166 

172 1F!AbS! Y2TEST) .LT.l .e--200)GO TO 173 
1F!AbS(Y(2)/Y2TEST)-1 .E-04 ) 1 73. 1 66. 1 66 

173 1F!ABS!Y3TEST).LT.1 .F-2001GO TO 174 

IF ! AbS ! Y ! 3 1/Y3TEST )-l .E-04 1174.166.166 

174 if!AbS!Y 4TEST1 .LT.l .P-2001GO TO 175 
1 F!AbS!Y!41/Y4TEST1-1 .E-04 1 175.166.166 

175 Y ! 1 1 =Y ! 1 1+Yl test 
Y !2 1 =Y ( 2 1+Y2TEST 
Y!31=Y!3)+Y3TEST 
Y!41=Y!4)+Y4TEST 


EMI=2.0 

EMJ=2.0 

1F<M1.EQ.01EM1 = 1 .0 
IF!MJ.EQ.01EMJ=1 .0 
MPLUE=M J+M I 

SORTXI =SQRT!XMNP(MIP,NI )**2-Ml**2 1 
SQRTxJ=SORT!XMNP (MJP ,NJ 1**2-MJ**2 1 

IF!TmI 1S0RTXI=1 .0 
IF!TMJ1S0RTXJ=1 .0 
MIP=Ml +1 

mup^mu^ 1 

CY=-sQRT (EM I *EMJ 1 / ! 1 20.*PI*SORTX I *S0RTX J 1 
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8899 


8891 

8801 

8802 

8803 

8804 
880?5 

8806 

8810 

3892 

97 


09999 
1 000 


1 002 


1 003 


1 004 


COEFF»CMPLX ( 0 , 0 . C Y ) 

CONTINUE 

YTE«COEFF*CMPLX<Y(l ) ,Y(3) ) 

YTM»COEFF*CMPUX(Y(2) «Y(4) ) 

IFIR.EO. 0,0)8891 « 889? 

1F(TM1 ♦AND, TMJ) 8801 *(3902 
VB*0,0 

UB=-(2,/EMl )*COS (FLOAT (MI )*PHIJPP) 

GO TO 8810 
IF(TmI >8803.8804 
VB=0,0 

UB=(FMI-1 ,0)*STN(FL0aT(MI )*PHIJPP) 

60 TO 8810 

IF (TmJ >8805. 8806 

V8=0,0 

UB*(FM1-1 ,0)*SIN(FLOAT(M1 )*PHIJPP) 

GO TO 8810 

VB=0»/EM1 )»COS (FLOAT (MI )*PHI JPP ) 

UB*-(EMI-1 ,0)*COS (FLOAT (Ml )*PH1 JPP) 

yte=yte*vb 

YTM3YTM#UB 

YC( I I . JJ)=YTE+YTM 

WR I TE ( 6 » 97 ) I I ♦ J J ♦ YC ( I I ♦ J J ) . PRMT ( 2 ) , 1 HOLE * 1 MODE . JHOLE . JMODE 
FORMAT! 1X*Y(*I2*.»12*) = 4E12.5* + J ( *E 1 2 , 5* ) *3X*PRMT ( 2 ) =*F6. 3 . 3X* I 
H0LEs*12.2X»lM0DE=*l2*3X»JH0LE=* I2.2X4JM0DE=*I2/ ) 

PH 1 JPP=PH 1 JPP* 1 80 , /P I 

RsR/OlMEN 

GO TO 1000 

YC ( I 1 . JJ)=CMPLX(0,0*0.0) 


continue 

DO 1001 IHOLE=l *NUMH0LE 
DO 1001 JHOLE=l .NUMHOLE 
DO 1001 IMODE=l .NUMMODE 
DO InOl JMODEal .NUMMODE 
TMI=, FALSE, 

IF( ImODE,GT,NUMTE)TMI=,TRUE, 

I 1 = ( I HOLE- 1 ) *NUMM0DE+ 1 MODE 
JJ=( JHOLE- 1 )*NUMM0DE+ JMODE 
AKZEROl =A I J ( IHOLE )*KZER0*D1 MEN 
1 F(TmI ) 1003.1002 
NI=NI J( IMODE) 

MIP=MI J( IM00E)+1 

FACT«;QI = (XMNP(MIP.NI )/AKZER0I )**2 
FTSQ=ER-FACTSQI 

IF(AbS(FTSO),LT.1 .OE-200)FTSO=0,0 
IF (FTSO.GE.0,0 ) YMNZERO= ( 1 ,0.0,0 )*SORT ( 
IF (FtSQ,LT,0,0 ) YMNZFRC=-CON*SQRT (-FTSQ 
GO TO 1 OO 4 


FTSQ) 

) 


I DEM= 1 MODE-NUMTE 
NI =NNl J ( IDEM ) 

MIP=mMI J( IDEM)+1 

FACTfOI = (XMN(MIP.NI )/AKZEROI )**2 

ftsq=er-factsqi 

CONVfRT=FTSQ 

IF(ApS(FTSO),LT,1 ,0E-290 )C0NVERT= 1 ,0E-290 
IF (FTSQ,GE.0,0 )YMNZEPO= ( 1 ,0.0,0 ) /SORT ( CONVERT ) 
IF(FtSO,LT,0,0)YMNZEpC= (-1 , .0, )/ (CON*SQRT (-CONVERT ) ) 
CONTINUE 

A( I I ,JJ)=YC( I I .JJ) 

B( I I .JJ)=-YC( I I . JJ) 
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96 


lF(t I.EQ.JJ)A( 1 !.JJ)=YMNZER0/(120.*PI )+YC( ! I«JJ) 

IF(II.EQ.JJ)B(n.JJ)=YMNZERO/( 120 .*Pn-YC(n.JJ) 

YC n 1 • JJ ) ( I 1 « JJ ) 


YMNZpRO = YMN'ZERO/ ( 1 20**PI ) 

Tcr/Tf WR I TE ( 6 • Q8 ) I I ♦ YMNZERO « I HOLE « I MODE 

= iei3.5* .J,.E.2.5...3X..H0LE..I2.3X.1"0DE=. 


I 


12 ) 

1001 CONTINUE 

WRITf=;(6*94 ) 

WRITf(6«813 ) 

813 FORMAT! 10X*SCATTERING f^ATRlX*) 

CAUL CX I NV ( A , M , B , M , DfTERM , I P I VOT , I NOEX , MAX , I SC ALE ) 

DO 0O2 1 = 1 » M 

DO 8n2 J=1 *M 
B<1«J)=CMPLX(0,«0.) 

DO 8->0 K=1 «M 

B(I*J)=B(1«J)+YC(1«K>*A(K,J) 


830 CONTINUE 
803 CONTINUE 

WR1TF<6« 1 05 ) 


803 


DO 805 I = 1 « M 
DO 8o5 J=1 »M 
XDX=CABS(B( I «J) ) 
IF(XDX.LT.5.E-16)G0 TO 810 
1 SOLATE=20.*ALOG1 0 (XDX ) 

XXXXl = A I MAG ( B ( I « J ) ) 


XXXX9 = REAL (B ( I . J) ) 

PHASp= ( 180 . /PI )*ATAN?<XXXX 1 ,XXXX 3 ) 
WRITF( 6.803 ) 1 . J.B( I . J) . isolate. PHASE 
FORMAT! 1 X*S(*I 2 *. * 12 *) = *E 12 . 5 * +J(*E 12 . 5 *) 


*3X.F9.4* 


1 DEG.* ) 


DB*3X.F8.3* 


GO TO 805 

, 810 WR I Tf ( 6 ♦ 806 ) I • J ♦ P n * J ) 

• 806 FORMAT! 1 X*S(*I 2 *»*I 2 *) 

805 CONTINUE 


♦E12.5* +J(»E12.5*)*3X*BEL0W -300 DB*) 
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WRITf(6»95) 
format ( 1 H 


) 


GO TO 1 

601 WRITF<6»88) ^ ^ 

88 format ( 1 HO^MODE SUBSCRIPT’S ARE OUT 

900 WRITr(6*901 ) 

901 FORMAT! 1 HI *1 OX* lOHENn OF JOB ) 


OF RANGE OF 


XMNP ARRAY*) 


STOP 

602 WRITf!6*99) 

99 format ! 1 X*MOOE 


SUBSCRIPTS FOR TM MOOES OUT OF RANGE OF XMN ARRAY*) 


STOP 

END 
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c 


c 


c 

c 


c 


c 


c 


A 


A 

B 

C 

D 


A 

B 

C 

o 

E 

F 

G 


A 

B 


7 

14 

21 

22 

23 

20 


35 

42 

49 

56 


63 

70 

77 

Bl 


84 
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SUBROUTINE FINDC(BETA*Y»0ERY) 

dimension BESSEL (21 ) ♦DERIVCS ) •DERY( 1 ) tEXTRA (8*8) *HOLD(8 ) ♦ 
INDEX ( 1 00 ) ,PARM(5 ) • RUST (8 ) f SAVE (500 ) tY( 1 ) 

DIMENSIONS for COMMON VARIABLES 
D I mens ion V3 (50 ) * V31 (50 ) f V32 (50 ) • V33 (50 ) « 

V4 (50 ) » V4 1 (50 ) % V42 (50 ) ^ V43 (50 ) t 
W4 (50 ) ^ W4 1 (50 ) » W42 (50 ) f W43 (50 ) « 

W3(50)«W31 (50 ) ^W32 (50 ) • W33 (50 ) t 
Z(50) •ZN(50) fXMNP(8*3) ♦XMN(0»3) 

common - dimensioned variables 

common V3*V31 1 V32*V33fW3*W31 • W32 • W33 • Z ♦ XMNP t XMN • 

V4»V41 fV42«V47» W4*W41 tWAa^WAS^ZN^ 

COMMON - undimensioned VARIABLES 
BSQfCCA •CCBfD?#KINDfL3tL4^M0STfTMI ,TMJ* 

NEW# NP3 • NP4 , RkrERR # RK3 * RK4 , TERMA , TERMB , TERMC ♦ TERMO # 

VI# V2# V3XI3#V4XI4# 

W1 fWlSQ# W2#W2SQ#W3X!3#W4XI4#XI l#Xl2#Xl3#Xl4f 
MI #MJfNl #NJ#MIP#MJP,AKZEROI # AKZEROJ t RKZERO # FACTOR I tFACTORJ* 
FACTSQI ♦FACTSOJfCOSPf C0SM#S INP#SINM,C0SPHI JfPHIJPP 
logical TMIfTMJ 

complex cargo # CCA f CCB t CCON # COS I NE # 

FlBOfFlPBOtFlPXll fFlXIl #F2PX I 2 #F2X I 2 • 
G1B0,G1PB0#G1PXI1 #G1XI 1 # G2PX I 2 t G2X I 2 • <1 ♦K2#SINE 
external LAYER3#LAYEP4 
IF(NFW )7# 105#7 

FIND OUT IF THIS BETA IS IN SAVE TABLE 
IF (BfT A -SAVE (LAST ) ) 1 4#84#35 

beta is less than last table value USED 
LAST=LAST-1 

IF<BFTA-SAVE (LAST ) )2l #84#28 
IF(LAST-1 >22*22# 14 
WRITf(6*23) 

format ( 10X#8HERR0R 23) 

GO TO 900 
NEXTrLAST 
LAST=LAST+1 
GO TO 63 

BETA IS GREATER THAN LAST TABLE VALUE USED 
IF (MOST-LAST >42*42 #40 
NEED=1 
GO To 109 
LAST=LAST+l 

IF(BfTA-SAVE (LAST > >56 •84*35 
NEXT=LAST-1 

at this POINT WE KNOW THAT BETA LIES BETWEEN 
SAVE(NEXT) AND SAVE(LAST) 

IF (ArS ( (BETA-SAVE (NEXT ) )/BETA )-l *E-6 >70 # 70 # 77 

LAST=NEXT 

GO TO 84 

IF (ArS ( (BETA-SAVE (LAST > >/BETA)-l *E-6 >84 # 84 ♦ 8 1 

NFEDrO 

GO To 109 

get integrand values from save table 

NOW= INDEX(LAST > 

DO 9l I «1 #4 
DERYd >=SAVE(NOW> 

now=now^-i 

CONTINUE 

RETURN 

beta is zero 
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105 IFlBE^TA.NE.O.O )GO TO 109 
DO 107 1=1*4 

DFPYf I 1=0, 

107 CONTINUE 
GO TO 403 

c calculate integrands 

RUST ARRAY DEFINED 
RUSTd )=R(XI ) 

C RUST(a)=S(Xl ) 

C RUST(3)=T(XI ) 

C RUST(4)=U(XI) 

C PUST(5)=R-(XI) 

C RUST(6) = S-(X1 ) 

C DUST(7)=T-(XI ) 

C DUST(fl)=U-(XI ) 

109 BS0=RETA#BETA 

ROOT = SORT(ABS(BSO-l • ) ) 

GO TO( 154 . 133* 154 » 1 1 ? ) *KIND 
C 

C case 4 - X14=XI3=XI2 

C 

LI 2 RUSTd )=1 . 

RtJST(2)=0. 

RUST(3)=1 . 

RUST(4 1=0, 

1=-(B=-TA-1 , ) 1 19,1 19. 1?6 
1 19 RUST (5 1=0, 

RUST (6 ) =-ROOT 
RUST (7) =0, 

RUST{8)=-ROOT 
GO TO 273 
1 26 RUST ( 5 ) =-ROOT 
RUST (6 1=0. 

RUST(7) =-ROOT 
RUST(8)=0, 

GO TO 273 
C 

C CASE 2 - XI4=XI3 

C 

133 ASSIGN 135 TO JAIL 
PARMd )=XI3 
PARM(2)=XI2 
PARM (3 )=-«5 
135 RUST( 1 ) = 1 • 

RUST (2 )=0, 

RUST<3)=1 • 

RUST (4 )=0, 

IF(BRTA-1 , )140, 140, 147 
140 RUST(5)=0, 

RUST(6)=-ROOT 
RUST ( 7 1 =-W3X I 3*ROOT 
RUST ( 8 ) =- V3X I 3*ROOT 
GO TO 249 
147 RUST(5)=-ROOT 
RUST(6)=0, 

RUST ( 7 ) =-V3X I 3*ROOT 
RUST(8)= W3XI3*ROOT 
GO TO 249 
C 
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C CASE 1 OR 3 - SET INITIAL VALUES FOR INTEGRATION 

C PROM XI4 TO XI3 

C 

154 PARM<1 )aXI4 
PARM(2)*XI3 
PARM(3)»-.5 
161 PARM(4)«RK4 
PARM(5)*0, 

DO 168 Iel»8 
OERIVI I )=.125 
168 CONTINUE 
RUSTll )*I , 

RUST(2)*0. 

RUST ( 3 ) a 1 , 

RUST(4 )=0. 

Ic-(BfTA-1 . )175,175, 1B2 
175 RUST(5)=0. 

RUST(6)=-ROOT 
RUST ( 7 ) =-W4X I4*R00T 
RUST ( 8 ) =-V4X I 44R00T 
GO TO 189 
182 RUST(5)=-ROOT 
RUST(6)=0, 

RUST ( 7 ) a-V4X I 4*R00T 
RUST (8) a W4XI4*R00T 
C INTEGRATE FROM XI4 TO XI 3 

189 CALL LRKS2<PARMtRUST,0ERIV*8«LAYER4«EXTRA ) 

IF(PaRM(5) )900«203. 196 

196 RK4=RK4*10. 

WRITF(6«197)RK4 

197 FORMAT! 1 OX. 39HERROR TOLERANCE FOR LAYER 4 INTEGRATION./ 

A 10X.21HHAS BEFN INCREASED TO. E13.5/) 

GO TO 161 

203 IF(KIN0-3)217.210.217 

case 3 - XI3a0 - SET LAYER 1 FUNCTIONS AND PROCEED 
TO FINAL CALCULATIONS 
210 F1B0=CMPLX(RUST( 1 ).RuST(2) ) 

GlB0aCMPLX(RUST(3).RuST(4) ) 

F1PB0=CMPLX(RUST(5) .rust ( 6 ) ) 

G1PB0=CMPLX(RUST(7) .RUST( 8 ) ) 

GO TO 300 

CASE 1 - ALL LAYPRS - SET INITIAL CONDITIONS 

217 PARM( 1 )axi3 
PARM(2 ) axi2 
PARM(3)a-,6 
DO 294 Iai,8 
HOLD ( I ) aRUST ( I 1 
224 CONTINUE 

ASS I ON 231 TO JAIL 
GO TO 245 
231 DO 278 1=1.6 

RUST! I laHOLD! I ) 

238 CONTINUE 

245 RUST ( 7 ) =TERMA*H0L0 ! 7 ) -TERMB*HOLD ! 8 ) 

RUST ! 8 ) =TERMA*H0LD!8 >+TERMB*HOLD !7 ) 

249 DO 252 1=1.8 
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OERlv< t )*.125 
252 CONTINUE 

PARM(4 )=RK3 
PARM(5)»0* 

C INTEGRATE FROM Xl3 TO XI2 

CALL LRKS2 ( PARM « RUST , DER I V . 8 f LAVER3 .EXTRA ) 

IF (PaRM (5 ) )900.273.2f9 

259 RK3=PK3*10. 

WRITf(6«260)RK3 

260 FORMAT! 1 OX. 39HERR0R TOLERANCE FOR LAYER 3 INTEGRATION./ 

A 10X.21HHAS BEP-N INCREASED TO. El 3,5/) 

GO TO JAIL. < 135.231 ) 

DO NOT PASS GO 
DO NOT COLLECT S200 
273 IF<XI4«EQ«0»0)GO TO 300 
DO 280 1=7.8 
HOLD ( I ) =RUST ( I ) 

280 CONTINUE 

RUST ( 7 ) »TERMC*HOLD ( 7 ) -TERMD*HOLD ( 8 ) 

RUST ( 8 ) =TERMC*HOLD ( 8 )+TERMD*HOLD ( 7 ) 

C algebra from XI2 TO XII 

F2XIP=CMPLX(RUST ( 1 ) .RUST (2 ) ) 

G2XIP=CMPLX(RUST<3) .RUST(4) ) 

F2PXI2 = CMPLX(RUST <5 ) .RUST (6 ) ) 

G2PXI2=CMPLX(RUST (7) .RUST (8 ) ) 

VMBSQ=V2-BSO 
I F { W2 ) 282 . 282 . 286 

282 IF(VM3SQ)284.283.285 

283 FlXIl=F2XI2-F2PXI2*D2 
G1 XI 1 =G2XI2-G2PXI2*D? 

FIPXI 1=F2PXI2 

G 1 PX 1 1 = CCA*G2PX 1 2 
GO TO 291 < 

284 P2=0.0 
02=SoRT (-VMBSQ) 

GO TO 287 

285 P2=SoRT(VMBSQ) 

02 = 0,0 

GO To 287 

286 ROOT=SQRT(VMBSO*VMBSO+W2SO) 

P2=S0RT(,5*(R00T+VMBS0) ) 

Q2=S0RT ( ,5* (ROOT-VMBEO ) ) 

287 K2=CMPLX(P2.-02 ) 

CARG0=D2*K2 

S I NE=CS I N (CARGO ) 

COS1nE=CCOS (CARGO ) 

CCON=SINE/K2 

FI XI 1 =F2X 1 2*C0S I NE-F2PX I 2*CC0N 
G1 XI 1 =G2X12*C0S1NE-G2PXI2*CC0N 
FIPXI 1 =K24F2XI2*S1NE+F2PXI2*C0SINE 
GlPXl 1 =K2»CCA*G2XI2*5INE+CCA*G2PX12*C0SINE 
C ALGEBRA FROM XII TO 0,0 

291 VMBS0=V1-BS0 
IF(W1 1292.292.296 

292 IF(VmBS0)294.293.295 

293 FlBO=FlXn-FlPXll*XI 1 
GiBOrGlXl 1-GlPXl 1 *XI 1 
F1PB0=F1PXI 1 
G1PB0=G1PXI 1 

GO TO 300 
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294 P!=0.0 

01 =S0RT (-VMBSO) 

GO TO 297 

295 PI *SORT (VMBSO) 

01 *0,0 

GO TO 297 

296 R00T=S0RT<VMBS0*VMBS0+W1S0) 

Pl*S0RT(*5*(R00T + VMBc;0) ) 

0 1 =SORT ( . 5* ( ROOT- VMB«^Q ) ) 

297 K1 =CmPLX(P1 ,-01 ) 

CARGO=XIl*Kl 

S I NE=CS I N ( CARGO ) 

COSInE=CCOS(CARGO ) 

CCON=SINE/Kl 

F 1 BOstF 1 X 1 1 *COS 1 NE-F 1 PX n *CC0N 
G1B0=G1XI1*C0SINE-G1PX1 1*CC0N 
CC0N=K1*S1NE 

F1PB0*F1XI l*CCON+FlPxI l*COSINE 
G1PB0=G1X1 1#CC0N+G1PX1 1*C0SINE 

C 

C FINAL integrand CALCULATIONS 

C 

C GET BESSEL FUNCTIONS NEFDEO 

BESSEL (1 )*JO(ARG) 

BESSEL (2 )=J1 (ARG) 

C 

300 CONTINUE 

IF (RkZERO.EO.O.O IGO TO 301 

r^PLUS*MJ+MI 

MMINUS*MJ-M1 

IF (MMINUS,LT,0)MMINUS=M1-MJ 

arg*rkzero*beta 

MMP=mPLUS+1 

CALL BESJSf ARG, BESSEL «MMP) 

BESP=BESSEL ( MMP ) 

MMP*MMINUS+1 

CALL BESJSI ARG, BESSEL •MMP) 

BESMsBESSEL ( MMP ) 

IF (MI »GT,MJ )BESM= (-1 ,0 ) *MM I NUS*BESM 
IF(TmI ,AND.TMJ)31 1 ,31 2 

311 VBETa=0.0 

UBETA*- (-1 *0)**MJ*(BESP*C0SP+(-l . 0 ) **M I *BESM*COSM ) 
GO TO 315 

312 IF(Tm1,0R,TMJ)313,314 

313 VBET4»0,C 

UBETA* (-1 ,0 )**MJ* (BESP*SINP- (-1 « 0 ) **M I *BESM*S I NM ) 
GO TO 315 

314 VBETA* ( -1 «0 (BESP*COSP+ (-1 «0 ) **M I *BESM»C0SM ) 
UBETA* ( -1 ,0 )»*MJ* (BESP*COSP- (-1 ,0 1 »*MI *8ESM*C0SM ) 

315 CONTINUE 
GO TO 302 

301 UBETA*1«0 
VBETA*1«0 

302 1F(TmI ) 1310, 1320 

1310 MMP=MI+1 
MMPP=MMP+1 

IF (ArS(FACTSOI-BSO)-I .OE-06) 1311,1312,1312 

1311 ARG*XMN(M1P,N1) 

CALL BESJS( ARG, BESSEL tMMPP) 

X IBEtA=BETA*BESSEL ( MmPP ) / (2 • *ARG ) 
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Y I beta =0.0 

GO TO 350 

1312 ARG=AKZEROI*BETA 

CALL BESJS(ARG.BESSEL«MMP) 

X IBETA=BETA*BESSEL ( MWP ) / (FACTSQ I -BSQ ) 
yibeta=o.o 

GO TO 350 

1320 MMP=MI+1 

IF< ArS(FACTSOI-BSQ)-! .OE-06 ) 1 32 1 , 1 322.1322 

1321 ARG = XMNP(MIP.N 1 ) 

CALL BESJS( ARG, BESSEL. MMP) 

YIBETA = BESSEL (MMP )♦ ( aRG*ARG-MI I ) / ( 2 . ^FACTOR I 1 
GO TO 1323 

1322 MMPP=MMP+1 
ARG=AKZER0I»BETA 

CALL BESJS (ARG. BESSEL •MMPP) 

Y I bet A = XMNP ( M I P . N 1 ) *FACT0R I * ( M I *BESSEL ( MMP ) /ARG-BESSEL ( MMPP ) ) / ( FAC 
1 TSOI -BSO) 

1323 XIBEtA=MI*BESSEL(MMP)/BETA 
350 IF (TmJ) 1330. 1340 

1330 MMP=MJ+1 
MMPP=MMP+1 

IF(ArS (FACTSQ J-BSO) -1 .OE-06 11331 .1 332. 1 332 

1331 ARG=XMN(MJP.NJ) 

CALL BESJS (ARG. BESSEL ♦MMPP) 

X JBEtA=BETA*BESSEL ( MmPP ) / ( 2. *ARG ) 

YJBETA=0.0 
GO To 360 

1332 ARG=AKZER0J*BETA 

CALL BESJS( ARG. BESSEL. MMP) 

XJBEtA=BETA*BESSEL (MMP 1/ (FACTSQJ-8SQ 1 

YJBETA=0.0 

GO TO 360 

1340 MMP=mJ+1 

IF(ArS(FACTSQJ-BSQ)-1 .OE-06) 1341 .1342.1342 

1341 ARG=XMNP(MJP.NJ) 

CALL BESJS(ARG.BESSEl«MMP) 

Y JBEtA=BESSEL ( MMP ) * ( aRG* ARG-MJ*M J ) / ( 2 . *FACT0R J 1 
GO TO 1343 

1342 MMPP=MMP+1 

arg=akzeroj*beta 

CALL BESJS( ARG. BESSEL 'MMPP) 

Y JBEtA =XMNP ( MJP « N J ) *FACT0R J* ( MJ*BES SEL ( MMP ) /ARG-BESSEL ( MMPP ) ) / ( FAC 
1 TSO J-BSO) 

1343 XJBETA=MJ*BESSEL (MMP )/BETA 
360 CONTINUE 

IF(TMl .OR.TMJ)VBETA=O.0 
FACTrM=X1BETA*XJRETA*UBETA*BETA 
FACTrE=-YIBETA*YJBETa*VBETA*BETA 
IF(XI4.EO.O,0)GO TO 3015 
C SEPARATE REAL ANn IMAGINARY PARTS 

NTYPf=0 

3010 IF(ArS(REAL(F1B0) ).Gt.1.E120)G0 to 3011 
IF(ArS( AIMAG (F lBO ) ) .OT.I .E120 )G0 TO 3011 
IF(ArS(REAL(F1PB0) ) .GT.1 .E120)G0 TO 3011 
IF(ArS(AIMAG(F 1PB0) i.GT.l ,E120)G0 TO 3011 
NTYPf=0 

GO TO 3012 

3011 NTYPF=NTYPE+ 1 
if(ntype.gt.5)GO to poo 
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FlBOaFlBO*( 1 .E-lSOtO, ) 

F1PBO=F1PBO*(1 ,E-120«0. ) 

GO TO 3010 

3012 IE(AbS(REAL(G1B0) ),GT*1«E120)G0 TO 3013 
IF<ArS(A1MAG(G1B0))«GT.1«E120)G0 to 3013 
IF(AbS{REAL(G 1PB0 > ) .OT.I *E120)G0 TO 3013 
IF(AgS(AlMAG(GlPBO) ),GT.1.E120)GO TO 3013 
NTYPpsO 
GO TO 3014 

3013 NTYPF»NTYPE+1 
if(ntype.gt.5)go to poo 
GlBO=GlBO»n .E-120.0, ) 

G1P80=G1PBO*(1 ,E-120,0. ) 

GO TO 3012 

3014 CONTINUE 
CCONaFlPBO/FlBO 
OERY { 1 ) *FACTRE*REAL < rCON ) 

OERY ( 3 ) »FACTRE*A 1 MAG (CCON ) 

CCON=CCB*Gl BO/Gl PBO 
OERY ( 2 > =FACTRM*REAL ( CCON ) 

OERY (4 ) *FACTRM*A 1 MAG (CCON ) 
IF(X|4.NE.0,0)G0 TO 3016 

3015 OERY( 1 )»FACTRE*RUST(5) 

OERY < 3 ) «FACTRE*RUST ( 6 ) 

0ERY(2 )*0»0 
OERY (4 ) =-FACTRM/R00T 
1F(BfTA.LE«1 *0)G0 TO 3016 
0ERY(2)*-0ERY<4) 

OERY (4 ) *0.0 

3016 CONTINUE 

PAVE INTEGRANOS 

IF (NFW) 449*403 »400 
400 IF(NFE0)414«414.407 
403 NEW=1 

C MOVE TO TOP OF BfTA SAVl 

407 M0ST=M0ST+1 
LASTaMOST 
GO TO 428 

C make space in MlnOLE OF 

414 MOSTxMOST+1 
NA=LAST+1 
MOVEaMOST+NA 
00 421 N=NA*MOST 
MxMOvE-N 

SAVE(M )»SAVEtM-l ) 

INOEXIM )*IN0EX(M-1 ) 

421 CONTINUE 

C SAVE BETA ANO POINTER 

428 SAVE (LAST )=BETA 

INDEX (LAST) = 101 +44 (MOST-1 ) 

NOW* INDEX (LAST) 

C PAVE INTEGRANOS 

00 435 1=1.4 
SAVE(NOW)aOERY( I ) 

N0W=N0W+1 
435 CONTINUE 

C CHECK FOR TABLE FULL 

IF ( 1 OO-MOST )442.442 .800 


TABLE 


BETA SAVE TABLE 
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442 NEW=-1 

GO TO 800 

BETA SAVE TABLE tS FULL 
449 KEEP= I NOEX ( 1 ) 

I F ( NPEO ) 456 » 456 ♦ 470 

PUSH DOWN SAVE TABLE FROM SAVE (NEXT) 

456 LIMIT=NEXT-1 

DO 4ft3 1=1 tLIMlT 
SAVE( 1 )*SAVE( I+l ) 
index! I )*INOEX< I+ l ) 

463 CONTINUE 
LASTsNEXT 
GO TO 484 

PUSH DOWN ENTIRE BETA SAVE TABLE 
470 DO 477 1=1 «99 

SAVE! I )*SAVE! 1+1 ) 

INDEX! 1 )=1NDEX! 1+1 ) 

477 CONTINUE 
LASTsMOST 

5AVE BETA AND POINTER 
484 SAVE (LAST )=BETA 
INDEX (LAST)=KEFP 

SAVE integrands 
DO 491 1=1*4 

SAVE(KEEP)=DERY( I ) 

KEEP=KEEP+1 
491 CONTINUE 
800 RETURN 

900 WR1TP(6«901 ) 

IF !NTYPE»NE«0 ) WRI TE ! 6* 904 )F1B0*F1PB0«G1B0 » G1 PBO 
904 format < 1 X*F 1B0 = *2E13«5/1X*F1 PB0 = *2E 13*5/1 X*G1 B0=*2E13*5/1 X*G 1 PBO 
12E13,5/). 

901 FORMAT! 1 HI « 1 Ox* lOHENO OF JOB ) 

STOP 

WRITp!6*903) 

format ( IX+MPLUS EXCEEDES 

stop 


902 

903 


DIMENSION OF ARRAY BESSEL*) 
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SUBROUT I NE SPLRED (N « pPSLN , X . Y * DEL Y « S2 t S3 ) 

DIMENSION X( n .Y( 1 ) .nELY( 1 ) ♦ S2 ( 1 ) t S3 (1 ) , 

A H ( 1 00 ) . H2 ( 1 00 ) «B ( 1 00 ) • DELSQY ( 1 00 ) ♦ C ( 1 00 ) 

N1 *N-1 
DO 7 1=1 *N1 

H( I )=x( i+i )-x( n 

7 DELY( I ) = (Y( n-Y( 1+1 ) )/H( 1 ) 

DO 14 l»2tNl 

H2( 1 )=H( I-l )+H< 1 ) 

B( 1 )s«5*H( 1-1 )/H2 ( 1 ) 

DEUSOYl 1 ) = (DELY( I )-DFUY( 1-1 ) )/H2 ( 1 ) 

S2 ( 1 )=DELSQY( 1 )+OELSOY( 1 ) 

14 C{I)=S2( 1 )+DELSQY(l ) 

S2(l 1=0. 

S2 (N)=0. 

0MEGA=1 .071797 
2! ETA=0. 

DO 3F 1=2, N1 

W= (C ( 1 )-B ( 1 )*S2( I-l )-( .5-B( I ) )*S2 ( I + l )-S2 ( I ) )*OMEGA 
lFtARS(W)-ETA)35,35»28 
28 ETA=aBS(W) 

35 S2 ( 1 )=S2( 1 )+W 

IF (EtA-EPSLN 142,21 ,2l 
42 DO 49 1=1 ,N1 

49 S3 ( I )= (S2 ( I+l 1-S2 ( I ) )/H( I ) 

RETUPN 

END 
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SUBROUT I NE SPL02 (N*M*T«X«Y«Z» DEL Y »DELZ * S2 *T3 « S3* T3 « SS*TT ♦ SSI * TT 1 ) 
DIMENSION X ( 1 ) «Y < 1 ) * DELY ( 1 ) * DELZ (1 )<S2(1 )«T2<1 )«S3(1 )«T3<1 )*Z(1 ) 
DATA SIXTH/. 166666666666667/ 

7 I=M 

IF (M-l )77*21 . 1 4 
14 1F(M-N )21 .21 *77 
21 IF(T-X( 1 ) 163. 28*35 
28 1 = 1 

GO TO 105 

35 IF(T-X(N) >42*91 *73 
42 1F(T-X( 1 ) >56*105*49 
49 1=1+1 

IF<T-X( I > >98* 105*49 
56 1=1-1 

IF(T-X( I ) >56*105*105 

63 IF<T-X( 1 >+l .E-6)65. 64*64 

64 T*X< 1 > 

* GO TO 28 

65 WR1TF<6*70) T 

70 FORMAT! 1 OX* IOHARGUMEnT = * E 1 4 . 6. 22HOUT OF RANGE IN SPLD2 ) 
WRITf(6*71 > (X(U)*L=1 *N) 

71 FORMAT(/10X*24HRANGE OF ARGUMENT VALUES/ ( E20.6 ) ) 

STOP 

73 1F(T-X(N)-1 .£-6)75*76*65 
*75 T=X(N> 

GO TO 91 
77 WRITF(6*84> 

84 FORMAT! 10X.23HM OUT OF RANGE IN SPL02 ) 

STOP 
91 I=N 
98 1=1-1 
105 HT1=T-X! I > 

HT2 = T-X ! 1 + 1 > 

PR0D=HT1*HT2 

SS2 = S2 ! I >+HTl*S3 ( I > 

TT2 = T2! I >+HTl*T3! 1 ) 

DELSqS= !S 2! 1 >+S2( 1+1 )+5S2)+SIXTH 
DELSQT= !T2! 1 ) + T2 ( 1 + 1 > + TT2 )*SIXTH 
SS=Y! I >+HTl*DELY! 1 > +PROD+DELSQS 
TT=Z ! 1 >+HT1*DELZ ( I ) +pROD*DELSQT 
H12=HT1 +HT2 
PRCOn=PROD*S 1 XTH 

SS1=DELY! I >+Hl 2+DELSQS+PRC0N+S3 ( 1 ) 

TT1=0ELZ! 1 )+H12*DELS0T+PRC0N*T3 ( 1 ) 

M= I 

return 

END 
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c 


c 


SUBROUTINE LAYERS (XI «RUST. DERI V ) 

DIMENSION RUST{ I ) »DERI V( 1 ) 

DIMENSIONS for COMMON VARIABLES 
dimension V3(50 ) ♦ V31 (50 ) » V32 (50 ) ♦ V33 (50 ) « 

A V4 (50 ) ♦VAl (50) »V42 ( 50 ) « V43 ( 50 ) ♦ 

B W4(50)»W41 (50).W42(50)«W43(50). 

C W3(50 ) ♦ W31 (50 ) *W32 ( 50 ) . W33 ( 50 ) ♦ 

D 2(50) .ZN(50)»XMNP(8*3) «XMN(8«3) 

COMMON - DIMENSIONED VARIABLES 
common V3«V31 ♦V32»V33»W3»V(31 « W32 t W33 ♦ 2 ♦ XMNP ♦ XMN * 

A V4»V41 »V42*V43»W4«W41 »W42*W43«2N« 

COMMON - UNoIMENSIONED VARIABLES 
B BSQ.CCA»CCB.D2»KIND»L3*L4«M0ST«TM1 «TMJ. 

C NEW . NP3 ♦ NP4 » RXERR . RK3 « R<4 ♦ TERMA « TERMB « TERMC « TERMD « 

D Vl» V2» V3X13«V4X14. 

E W1 «W1S0.W2»w2SQ»W3XI3«W4X14.X1 1«XI2.XI3«XI4, 

F MI .MJ»NI ,NJ.MIP«MJP,AK2ER0I .AKIERO J »RK2ER0 . FACTOR 1 .FACTORJ* 

G FACtSOI tFACTSOJ»COSPtCOSM»S I NP « S I NM * COSPH IJ ♦ PH 1 JPP 
logical TMI«TMJ 
complex ccAtCce 

CALL SPLD2(NP3.L3.XI ,2.V3.W3«V31 .W31 . V32 ♦ W32 « V33 ♦ W33 . V3X I ♦ 

A W3XI «V3PXI .W3PXI ) 

DERIVd )*RUST(5) 

DERIV(2)=RUST(6) 

DERIV(3)»RUST(7) 

DERIV(A )*RUST(8 ) 

BS0MV®BS0-V3XI 

DER I V ( 5 ) «BSQMV*RUST ( 1 ) -W3X I *RUST ( 2 ) 

0ERIV(6 )*BSOMV*RUST (? )+W3XI *RUST ( 1 ) 

DEN0M=V3XI*V3XI+W3XI*W3XI 

FIRST=V3XI/DEN0M 

second® V3PX I *RUST ( 7 ) +W3PX I 4RUST ( 8 ) 

TH IR0*W3X I /DENOM 

F0URtH=V3PXI *RUST (8 )-W3PX I *RUST ( 7 ) 

DER I V ( 7 ) =F I RST*SECONd-TH I RD*FOURTH+BSQMV»RUST ( 3 ) -W3X I *RUST ( 4 ) 
DERIV(8 )®FIRST*F0URTH+THIRD*SEC0ND+BSQMV*RUST (4 )+W3Xl*RUST (3 ) 

return 

END 
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SUBROUTINE LayER4 (X I ,RUST . OER I V ) 
dimension rust ( 1 ) «DERI V( 1 ) 

C DIMENSIONS FOR COMMON VARIABLES 

DIMENSION V3 (50 ) « V3l (50 ) «V32 (50 ) ♦ V33 ( 50 ) . 

A V4(5C)«V41 (50 ) «V42 (50 ) « V43 (50 ) « 

B W4(50)«W41 (50) «W42 (50 ) ,W43 (50 ) t 

C W3(50)«W31 (50) » W32 (50 ) »W33 (50 ) « 

D Z ( 50 ) . 2N ( 50 ) « XMNP (8»3).XMN(8»3) 

C COMMON - DIMENSIONED VARIABLES 

common V3.V31 .V32«V33«W3»W31 « W32 . W33 »2 . XMNP * XMN. 

A V4«V41 .V42«V43»W4»W4I «W42*W43*ZN» 

C COMMON - undimensioned VARIABLES 

B BSQ.CCA,CC3.D?»KIND»L3.L4.M0ST«TMI .TMJ, 

C NEW . NP3 * NP4 . RkERR.RK3.RK4 . TERMA , TERMS « TERMC . TERMD . 

D VI. V2. V3XI3.V4XI4. 

E W1 .W1 SC3.W2.W2SQ.W3X I 3. WAX I 4. XI 1.XI2.XI3.XI4. 

F MI .MJ.NI ,NJ.MIP.MJP,AKZEROI .AKZEROJ .RKZERO . FACTOR I .FACTORJ. 

G fact SO I .FACTSOJ.COSP.COSM.S INP. S INM .COSPHI J.PHI JPP 

logical tmi .tmj 
complex cca.ccb 

CALL SPLD2(NP4.L4.XI , ZN . V4 . W4 . V4 1 »W41 . V42 . W42 . V43 . W43. VAX I . 

A WAX I . V4PXI .W4PX1 ) 

DERIV( 1 )=RUST(5) 

DERIV(2 )=RUST(6) 

DERIV(3)=RUST(7) 

DERIV(4 )=RUST (8 ) 

BS0MV*BSQ-V4XI 

DER I V ( 5 ) =BSOMV*RUST ( 1 ) -WAX I *RUST ( 2 ) 

DER I V ( 6 ) =BSQMV*RUST ( 2 ) +W4X I *RUST ( 1 ) 

DENOM = VAX I *V4X I +W4X I *W4X I 
FIRST=V4X I/DENOM 

second* V4PX I *RUST (7 ) +W4PX I *RUST ( 8 ) 

THIRn=W4XI /DENOM 

FOURTH* V4PX I *RUST (8 ) -W4PX I *RUST ( 7 ) 

DER I V ( 7 ) =F I RST*SEC0ND-TH 1 RD*FOURTH+BSQM V*RUST ( 3 ) -WAX I *RUST ( A ) 
DER I V ( 8 ) =F 1 RST*F0URTH+TH I RD*SECOND+BSQMV*RUST ( 4 ) +WAXI *RUST ( 3 ) 
RETURN 
END 
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SUBROUTINE BES JS ( XX . RJ .MT ) 

DIMENSION BJ(1)»B(130) 

C ROUTINE FINDS BESSEL J OF X FOR ORDERS ZERO THROUGH MT 

C AND loads them InTO BJ(1) THROUGH BJ(MT+1). 

X=ABS(XX) 

BJ ( M = 1 *0 
N * MT + 1 
IF (X.GE»80* ) I 0,20 

10 PI =3, 14 1592653589793 
DO 11 I «1 ,N 

11 BJ( I )=SQRT(2*0/<PI*X) )*COS<X-0»25*PI-0,5*P1*( I-l ) ) 

GO TO 220 

20 CONTINUE 

DO 5 I = 2»N 
5 BJ ( I ) = .0 

IF(X-15. )32,32,34 
32 NTEST = 20,+l0,*X-X»*2/3 
GO TO 36 

34 NTEST = 90.+X/a, 

36 IF(MT-NTEST)40, 38,38 
38 N = nTEST - 1 
GO TO 45 
40 N ® MT 
45 8PREV = «0 
N1 s N+1 
F = ?,/X 
0 » l.OE-6 

C COMPUTE STARTING VALUE OF M 

IF (X-5. >50,60,60 
50 MA = X + 6» 

GO TO 70 

60 MA s 1,4»X + 60,/X 

70 MB * N+IFIX(X)/4+2 

MZERO = MAX0(MA,MB> 

C SET upper limit OF M 

MMAX = NTEST 
DO igO M = MZERO, MMAX, 3 
FMl = l,0E-28 

FM = ,0 
ALPHA = .0 

IF (M- (M/2 )*2 ) 1 20, 1 1 0, 1 20 
1 1 0 JT = - 1 

GO TO 130 
120 JT = 1 

1 30 M2 = M-2 

DO 160 K = 1 ,M2 

MK = M-K 

B(MK) = F*FL0AT(M<)*FM1-FM 
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overflow test 

tF (B (MK )-l .OEeS ) 140»?20«2a0 
140 FM = FMl 

FMl = 8(MK) 

JT = -JT 
S = 1+JT 

160 alpha = ALPHA+B(MK)*c 
0(1) = f*fm1-FM 
alpha = ALPHA+B(1) 

BTEST = B (N1 ) 

BTEST = BTEST/ALPHA 

IF(AbS{BTEST-BPREV)-ABS(D*BTEST ) ia00«200* 190 
190 BPREV = BTEST 
200 00 2l0 1 = 1 »N1 

Bin BJU) = B(1)/ALPHA 
220 1 F (XX»LT.0.0 )G0 TO 2^0 
RFTUPN 

230 N=MT+1 

00 201 1=1«N 

NN= I - 1 

231 BJ( I >=8J( I )*(-! •0)**NN 
RETURN 

END 
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SUBROUTINE LRKSI ( PRMt • Y * DERY » ND I M «FCT « AUX ) 

dimension Y( 1 > «DERY( 1 ) « AUX(8«2) «A (4 ) ♦BIA) «C(4 ) «PRMT( 1 ) 

50 format (52H MORE THAN 15 BISECTIONS NEEDED IN LRKSI INTEGRATION) 

51 F0RMAT(47H initial INCREMENT IS ZERO ON LRKSI INTEGRATION ) 

52 F0RMATI54H INITIAL INCREMENT HAS WRONG SIGN IN LRKSI INTEGRATION) 
DO I lal »NDIM 

1 AUX(fl*I) = .06666666A6666667*DERY( I ) 

XsPRmT(I) 

XEND=PRMT (2 ) 

HsPRmT ( 3 ) 

CALL FCTIX.Y.DERY ) 

ERROR TEST 

I F ( H* ( XEND-X ) ) 38 t 37 t ^ 

preparations for runoe-kutta method 

2 AM )=.5 

A (2) = .292893218813452 

A(3) = 1.70710678118655 

A (4) = .166666666666667 

B(1 )=2. 

B(2)sl. 

B(3)=l . 

B(4)s2. 

CM )s.5 
C(2) = A(2) 

C(3) = A(3) 

C(4)s.5 

preparations of first runge-kutta step 

do 3 I =1 .NDIM 
AUX( 1 . I )=Y( 1 ) 

AUX(2. I )=0ERY{ I ) 

AUX(7. I )=0. 

3 AUX(6* I 1*0. 

IREC=0 

H=H+H 

IHLF=-1 

I STEP“0 
IEND=0 


START OF A RUNGE-KUTtA STEP 

4 IF( (X+H-XENn)*H)7.6«'=: 

5 H=XEND-X 

6 IEND=1 

7 ITEST = 0 

9 ISTEP=ISTEP+1 


start of innermost RuNGE-KUTTA loop 
J = 1 

10 AJ=A(J) 

BJ=B( J) 

CJ=C< J) 

DO 11 1=1. NDIM 

R1 =H»DERY ( I ) 

R2=AJ*(R1 -BJ*AUX(6. I ) ) 

Y( I ) =Y( I )+R2 
P2=R?+R2+R2 
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11 AUX(6» I )=AUX(6* I )+R2-CJ*Rl 
IF( J-4 1 12. 15* 15 

12 J=J+1 

lF(J-3) 13.14.13 

13 X»X+,5»H 

14 CALL FCTIX.Y.DERY ) 

GOTO 1 0 

C END OF INNERMOST RUNGE-KUTTA LOOP 

C 

C 

C TEST OF ACCURACY 

15 IF(ITEST)16.16.20 

IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 

16 DO 17 1 =1 .NDIM 

17 AUX(4. I )=Y( I ) 

I TEST” 1 

1 STEP” I STEP+ 1 STEP-2 

18 1HLF=1HLF+1 
X=X-H 
H=»5*H 

DO 19 1=1 .NDIM 
Y( I )=AUX( 1 . I ) 

DERY( 1 )=AUX(2. I ) 

19 AUX(6. I )=AUX(3. 1 ) 

GOTO 9 

IN CASE 1TEST=1 TESTING OF ACCURACY IS POSSIBLE 

20 IMOD=ISTEP/2 
I F ( I ^TEP- 1 MOD- I MOD ) 2 1 . 23 . 2 1 

21 CALL FCTfX.Y.DERY ) 

DO 2? 1=1 .NDIM 

AUX(«^. 1 )=Y( I ) 

22 AUX(7. I )=DERY( I ) 

GOTO 9 

computation of test value DELT 

23 DELT=0, 

DO 24 1 =1 .NDIM 
IF (Y( 1 ) 1 242. 241 . 742 

241 DELT = DELT + AUX(8.I> * ABS ( AUX ( 4 . 1 ) ) 

GO To 24 

242 DELT = DELT + AUX (8. I) * ABS ( ( AUX ( 4 . I ) - Y(I)) / Y(I)1 

24 CONTINUE 

IF(DFLT-PRMT(4) 128.28.25 

error IS TOO GREAT 

25 1F( lHLF-15126.36.36 

26 DO 27 I =1 .NDIM 

27 AUX (4. I )=AUX(5. 1 1 
I STEP” I STEP+ 1 STEP-4 
X=X-H 
IEND=0 
GOTO 18 

C result values are GOOD 

28 CALL FCT(X.Y.DERY) 

DO 29 1=1 .NDIM 
AUX( 1 . I )=Y( 1 1 

AUX (2. I )=DERY( I ) 
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AUX(?. 1 ) =AUX(6* 1 ) 

Y( I )sAUX(5» I ) 

29 OERY( I »*AUX(7. I ) 

30 DO 31 1 =1 .NDIM 

Y( I IsAUXd . I ) 

31 OERYf I )=AUX(2, I ) 

IREC=1HLF 
IF( IFN0)32«32.40 

INCRFMENT GETS DOUBLED 

32 IHLFsIHLF-1 
ISTED=1STFP/2 
H=H+h 

IF( 1HLF14, 33*33 

33 I MOD=l STEP/2 
TF (I STEP- I MOD- I MOO 14 ,34,4 

34 IF (DFLT-,02*PRMT (4 ) ) S5,35,4 

35 IHLF=IHLF-1 
ISTEp=ISTFP/2 
H = H+h 
GOTO 4 

c 

c 

c returns to calling program 

36 WRITf(6»50) 

PRMT (5 ) =1 ,0 
GO TO 40 

37 WRITP(6*51) 

GO TO 39 

38 WRITp(6.52) 

39 PRMT(5)=-1 .0 

40 return 
END 


Continued 
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SUBRoUT I NE LRKS2 ( PRMT ♦ Y « DER Y « ND I M ♦ FCT t AUX ) 

DIMENSION Y( 1 ) *DERY ( 1 ) « AUX (8 « 1 )«A(4)»B(4)«C(4) »PRMT ( 1 ) 

50 F0RMAT(52H more THAN 15 BISECTIONS NEEDED IN LRKS2 INTEGRATION) 

51 F0RMAT(47H initial INCREMENT IS ZERO ON LRKS2 INTEGRATION ) 

52 format (54H INITIAL INCREMENT HAS WRONG SIGN IN LRKS2 INTEGRATION) 

DO 1 I = 1 « ND I M 

1 AUX(8»I) = •0666666666666667*DERY( I ) ^ 

X=PRmT ( 1 ) 

XEND=PRMT (2 ) 

H=PRmT ( 3 ) 

CALL FCT(X.Y,DERY ) 

C 

c error test 

IF(H*(XEND-X) )38«37»2 
C 

C preparations for RUNGE-KUTTA METHOD 

2 A(1 )=.5 

A(2) = .292893218813452 

A(3) = 1.70710675118855 

A (4) = .166666666666667 

B(1 )=2. 

B(2)=l. 

B(3)=l . 

B(4)=2. 

C(1 )=.5 
C (2 ) = A (2 ) 

C(3) = A(3) 

C (4 ) = .5 
C 

c preparations of first runge-kutta step 

DO 3 I = 1 . ND I M 
AUX< 1 . I )=Y( I ) 

AUX(2. I )=OERY( I ) 

AUX(3« I )=0. 

3 AUX(6« I )=0. 

IREC=0 

H=H+H 

IHLF=-1 

ISTEP=0 

iend=o 

c 

c 

C start OF A RUNGE-KUTtA STEP 

4 IF( (X+H-XFN0)*H)7»6»5 

5 H=XEND-X 

6 IFNDsl 

7 ITESt = 0 

9 ISTEP=ISTEP+1 
C 
C 

c start of innermost Runge-kutta loop 
j=i 

10 AJ=A(J) 

BJ=B(J) 

CJ=C( J) 

DO 11 1*1.N01M 

R1 =H*DERY( 1 ) 

R2=AJ*<R1-BJ*AUX(6. I ) ) 

Y( I )=Y( I )+R2 
R2*R2+R2+R2 
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11 AUX<6* I >=AUX<6« I )+R2-CJ*Rl 
IF(J_4)12«15.15 

12 J»J+1 
IF(J-3)13«14«13 

13 X»X+,5»H 

14 CAUL FCT(X»Y.D<ERY ) 

GOTO 10 

C END OF innermost RUNrE-KUTTA LOOP 

C 

C 

C TEST OF ACCURACY 

15 IF ( ITEST 116. 16.20 

IN CASE ITESTsO THERF IS NO POSSIBILITY FOR TESTING OF ACCURACY 

16 OO 17 1=1 .NDIM 

17 AUX(4» 1 )=Y( 1) 

ITEST=1 

I STEP= I STEP+ 1 STEP-2 

18 1HLF=IHLF+1 
X«X-H 
H=«5»H 

DO 19 I=1.N01M 
Y( I )=AUXU • I ) 

DERY( 1 )«AUX{2. 1 ) 

19 AUX(6. I )=AUX(3. I ) 

GOTO 9 

IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 

20 IMOD=ISTEP/2 
I F ( I FTEP- I MOD- I MOD ) 2 1 . 23 . 2 1 

21 call FCTIX.Y.DERY ) 

DO 22 1=1 .NDIM 
AUX(.?. 1 )=Y(1 ) 

22 AUX(7. I >«DERY( 1 ) 

GOTO 9 

computation of test value DELT 

23 DELTsO, 

00 24 I = 1 . NO I M 
IF (Y(IM 242. 241, 2A2 

241 DELT = DELT + AUXIB.t) * ABS (AUX(4, 1)) 

GO TO 24 

242 BELT = DELT + AUXIB.I) * ABS < ( AUX ( 4 . I ) - Yd)) / Yd)) 

24 CONTINUE 

1F(DPLT-PRMT(4) )28.2fl,25 
ERROP IS TOO GREAT 

25 IFdHLF-15)26.36.36 

26 DO 27 1 =1 .NDIM 

27 AUX(4. I )=AUX(5, I ) 

1 step* I STEP+ 1 STEP-4 
X=X-H 
IEND=0 
GOTO 18 

result values are GOOD 

28 CALL FCTIX.Y.DERY ) 

DO 29 I =1 .NDIM 
AUX( 1 , I )=Y( 1 ) 

AUX (2. 1 )=DERY( I ) 
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AUX(T. I )=AUX(6« I ) 

Yd )=AUX (5. I ) 

39 DERY( I )=AUX(7, I ) 

30 DO 3l 1=1 »NDIM 
Y ( I )=AUX < 1 d ) 

31 OERY( 1 )=AUX(3« 1) 

IREC* IHLF 

IF( lFND)33t33.40 

INCRrMENT GETS DOUBLED 
33 1HLF=1HLF-1 
I STEPnl STEP/3 
HaH+H 

1F(1hLF)4«33.33 

33 IMOD=ISTEP/3 

IF( 1 step- 1 mod- I MOD) 4, 34, 4 

34 IF(D=-LT-,03*PRMT(4) )-^B,35,4 

35 IHLF=1HLF-1 
lSTEP=lSTFP/3 
H=H+H 

GOTO 4 


RETURNS TO CALLING PROGRAM 

36 WR1Tp(6«50) 

PRMT(5)=1 ,0 
GO TO 40 

37 WR1Tp(6«51 ) 

GO TO 39 

38 WR1TP(6«52) j 

39 PRMT(5)=-1 .0 I 

40 return I 

END / > 

s 
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SUBROUTINE CXINV ( A «N ,8 »M »DET » IPI V» INOX* MAX* I SCALE ) 

COMPLEX MATRIX INVERSION WITH SOLUTION OF LINEAR EQUATIONS 

CAVM » CABS(A(MAX) ) * CAVA = CABS(A(I*J)) 

CAOM » CABS(DETERM) * CAPV « CABS (PIVOT) 

complex A(MAX*N)* B(MAX*M)* swap* OET * PIV* PIVl* CO* Cl 
DIMENSION IPIV(N)* IN0X(MAX*2) 

CONSTANTS* initialization 

CO = (0, 0*0,0) 

Cl = (1, 0*0,0) 

I scale = 0 

RL ■ 10,0**100 
RS = 1.0/RL 
OET s Cl 
CAOM = 1,0 

00 20 J“1 *N 
20 IPIV(J) * 0 

DO 500 1=1 *N 

SEARCH FOR PIVOT FLEMENT 

CAVM = 0,0 
DO 105 J*1*N 

IF (IPIV(J) ,E0, 1) GO TO 105 

DO 100 K=1 *N 

IF (IPIV(K) - 1) 50*100*750 
50 CONTINUE 

CAVA = CABS(A(J*K)) 

IF (CAVM ,GE, CAVA) GO TO 100 

1 ROW * J ' 

ICOL * X 

CAVM = CAVA 
100 CONTINUE 
105 CONTINUE 

IF (CAVM ,EQ, 0,0) GO TO 720 
IPIV(ICOL) = IPIVnCOL) + 1 

Interchange rows to put pivot element on diagonal 

IF (IROW ,E0. ICOL) GO TO 230 
OET = -OET 
00 200 L=1 *N 
SWAP = A(IR0W*L) 

A(IROW*L) = A(ICOL*L) 

A(ICOL*L) = SWAP 
200 CONTINUE 

IF (M ,LE. 0) GO TO P30 
00 2?0 L=1 *M 
SWAP = B(IROW*L) 

B(IROW*L) = B(ICOL*L) 

B(ICOL*L) = SWAP 
220 CONTINUE 
230 continue 

lNOX( 1*1) = IROW 
INDX(I*2) = ICOL 
PIV = A (ICOL* ICOL) 
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CAPV s CABS(PIV) 

IF (CAPV «EO. 0.0) GO TO 720 

c 

c Scale determinant 

c 

PIVI = PIV 
CADM = CABSIDET) 

IF (CADM ,LT, RL ) r-0 TO 360 
OET = DET/RL 
CADM = CABSIDET) 

ISCAlE = ISCALE + 1 
I-F - IfrA€>M~- .t:-T.- RL V- GO- TO 290 
OET = OET/RL 
ISCAlE = ISCALE 1 
60 TO 290 
260 CONTINUE 

IF (CADM .GT. RS) GO TO 290 
DET = DET*RL 
CADM = CABSIDET) 

ISCALE = ISCALE - 1 
IF ICADM ,GT. RS) GO TO 290 
DET s DET-»RL 
ISCAlE = ISCALE - 1 
290 CONTINUE 

CAPV = CABSIPIVI ) 

IF ICAPV ,LT. RL) GO TO 320 
PIVI = PIVI/RL 
CAPV = CABSIPIVI ) 

ISCALE = ISCALE + 1 

IF IcAPV ,IlT. RL) GO TO 340 
PIVI = PIVI/RL 
ISCALE = ISCALE + 1 
GO Tn 340 I 

,320 CONTINUE | f >' 

IF IcAPV »GT, RS) GO TO '340 
PIVI = PIVI»RL 
CAPV = CABSIPIVI ) 

ISCAlE = ISCALE - 1 
IF IcAPV .GT. RS) GO TO 340 
PIVI = PIVI*RL 
ISCALE = ISCALE - I 
340 CONTINUE 

DET = DET ♦ PIVI 
C 

C DIVIDE PIVOT ROW p>Y PIVOT ELEMENT 

C 

A I I COL* I COL ) = Cl 

DO 3f0 L=1 »N 

350 AIICOL.L) = A I ICOL.L )/PIV 
IF IM .LE. O) GO TO -<80 
DO 370 L=1 .M 

370 B(ICOL.L) = B I ICOL.L )/PIV 
C 

C Rc’DUCE NON-PIVOT DOWS 

C 

380 CONTINUE 

DO 500 LI =1 .N 

IF ILI .EQ. ICOL) GO TO 500 
SWAP = A ILI « ICOL ) 

A ILI . ICOL ) = CO 
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n n n 


APPENDIX - Concluded 


DO 4O0 L»1 »N 

400 A(Ul.L) = A(U1.L» - AnCOL«L)*SWAP 
IF (M .LF» 0) GO TO fOO 
DO 4F0 L=1 »M 

450 B(Ll.L) = B(L1*L) - B ( I COL *L )*SWAP 
500 CONTINUE 

interchange columns 

DO 7O0 I=1»N 
L F N+l-I 

IF (INDX(L«1) .EO. InDX(L.2)>G0 TO 700 

I ROW = INOXIL. I ) 

I COL » INOX (L *2) 

DO 6P0 Ksl.N 
SWAP = A(K«IROW) 

AIK. I ROW) = AIK. I COL) 

AIK.ICOL) = SWAP 
690 CONTINUE 
700 CONTINUE 
GO TO 750 
720 OET = CO 
. I scale = 0 
750 return 
END 
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